GAMES103 - Intro to Physics-Based Animation

(Based on Unity, C#)

Huamin Wang (games103@style3D.com) Video; Lecture site

Lecture 1 Introduction

Graphics

Geometry

Mesh: Triangle mesh is the foundation of graphics
Vertices (nodes) + Elements (triangles, polygons, tetrahedra…)
Point Cloud: simple, can be raw data from surface scan.
But problems in mesh reconstruction, (re)-sampling, neighborhood search …
Grid: often acquired from volumetric scan, e.g., CT
Problems in memory cost (octree?), volumetric rendering

Rendering

Photorealistic Rendering (Ray Tracing)
Non-Photorealistic Rendering

Material Scan

Body Scan by a Light Stage

Animations

Character Animation
Physics-Based Animation

Physics-Based Animation Topics

Rigid Bodies [contact / fracture] - Mesh / *Particle in fracture (to avoid remeshing)
Cloth and Hair [clothes / hair] - Mesh / *Grid (to simplify contacts)
Soft Bodies [elastic / plastic] - Mesh
Fluids [smoke / drops and waves / splashes] - Mesh in drops and waves (RT) / Particle (RT) in smoke and splashes / Grid (universal)

Hybrid Methods and Coupling - problems

Lecture 2 Math Background

Vector, Matrix and Tensor Calculus

Vectors

Definitions

A geometric entity endowed with magnitude and direction

\begin{matrix} p = [\begin{matrix} p_{x} \\ p_{y} \\ p_{z} \end{matrix}] \in R^{3}; o = [\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}] \end{matrix}

Right-Hand System (OpenGL, Research, …)
Left-Hand System (Unity, DirectX, …) => screen space

Can also be stacked up to form a high-dim vector -> describe the state of an obj (not a geometric vector but a stacked vector)

Arithematic

Addition and Substraction

(commutative)

\begin{matrix} p \pm q = [\begin{matrix} p_{x} \pm q_{x} \\ p_{y} \pm q_{y} \\ p_{z} \pm q_{z} \end{matrix}]; p + q = q + p \end{matrix}

Example $\vb{p}(t) = \vb{p} + t\vb{v}$ $0 < t<1$ $0<t$ $t\in \R$ $t$ is an interpolant)

Vector Norm

Magnitude of a vector (length)

1-norm $||\vb{p}||_1=|p_x|+|p_y |+|p_z|$
Euclidean (2) norm $||\vb{p}||_2 = (p_x^2 + p_y^2 + p_z^2) ^{1/2}$
p-norm $||\vb{p}||_p = (p_x^p + p_y^p + p_z^p) ^{1/p}$
Inifinity norm $||\vb{p}||_\infty=\max(|p_x |,|p_x |,|p_x |)$

Usage:

Distance $\vb{p}$ $\vb{q}$ $||\vb{q}-\vb{p}||$
Unit Vector $||\vb{p}|| = 1$
Normalization $\bar {\vb{p}} = \vb{p} / ||\vb{p}||$ $||\bar {\vb{p}}|| = ||\vb{p}|| / ||\vb{p}|| = 1$

Dot Product

(inner product)

< p, q >\equiv p \cdot q = p_{x} q_{x} + p_{y} q_{y} + p_{z} q_{z} = p^{T} q = | | p | | | | q | | \cos θ

Properties:

$\vb{p}\cdot \vb{q} = \vb{q}\cdot \vb{p}$
$\vb{p}\cdot \vb{(q+r)} = \vb{p}\cdot \vb{q} + \vb{p}\cdot \vb{r}$
$\vb{p}\cdot \vb{p} = ||\vb{p}||^2_2$ , an alternative way to write norm
$\vb{p}\cdot \vb{q} = 0$ $\vb{p}, \vb{q} \neq 0$ $\cos \theta = 0$ => orthogonal

Example: Particle-Line Projection

$s = ||\vb{q} - \vb{o}||\ \cos \theta$ $\vb{s} = \vb{o} - s\bar{\vb{v}}$

\begin{aligned} s & = | | q - o | | \cos θ \\ = | | q - o | | | | v | | \cos θ / | | v | | \\ = (q - o)^{T} v / | | v | | \\ = (q - o)^{T} \bar{v} (normalized) \end{aligned}

Example: Plane Representation

$s$ $\vb{n}$ - normal)

\begin{matrix} s = (p - c)^{T} n {\begin{aligned} > 0 & Above the plane \\ = 0 & On the plane \\ < 0 & Below the plane \end{aligned} \end{matrix}

Example: Particle-Sphere Collision

\begin{aligned} | | p (t) - c | |^{2} & = r^{2} \\ (p - c + t v) \cdot (p - c + t v) & = 𝑟^{2} \\ (v \cdot v) t^{2} + 2 (p - c) \cdot v t + (p - c) \cdot (p - c) - r^{2} & = 0 \end{aligned}

$t$ $t$ )

Cross Product

\begin{matrix} r = p \times q = [\begin{matrix} p_{y} q_{z} - p_{z} q_{y} \\ p_{z} q_{x} - p_{x} q_{z} \\ p_{x} q_{y} - p_{y} q_{x} \end{matrix}] \end{matrix}

Properties

$\vb{r} \cdot \vb{p} = 0$ $\vb{r} \cdot \vb{q} = 0$ $||\vb{r}|| = ||\vb{p}|| \ ||\vb{q}|| \ \sin\theta$
$\vb{p}\times\vb{q} = - \vb{q}\times\vb{p}$
$\vb{p} \times(\vb{q}+ \vb{r}) = \vb{p}\times\vb{q} + \vb{p}\times\vb{r}$
$\vb{p}\times\vb{q} = \vb{0}$ $\vb{p}, \ \vb{q} \neq 0$ $\sin\theta = 0$ $\vb{p}$ $\vb{q}$ are parallel (direction can be opposite)

Example: Trinagle Normal and Area

$\vb{x}_{10} = \vb{x}_1 - \vb{x}_0$ $\vb{x}_{20} = \vb{x}_2 - \vb{x}_0$
$\vb{n} = (\vb{x}_{10} \times\vb{x}_{20}) / ||\vb{x}_{10} \times\vb{x}_{20}||$ (dep on the topological order)
$A = ||\vb{x}_{10}|| h /2 = ||\vb{x}_{10} \times\vb{x}_{20}||/2$

Example: Trianle Inside / Outside Test (Co-plane)

$\vb{x}_0 \vb{x_1}$ $(\vb{x}_0 -\vb{p})\times (\vb{x}_1 - \vb{p})\cdot \vb{n} > 0$ (Same normal as the main triangle)
$\vb{x}_0\vb{x}_1$ $(\vb{x}_0 -\vb{p})\times (\vb{x}_1 - \vb{p})\cdot \vb{n} < 0$

Example: Barycentric Coordinates

\begin{matrix} \frac{1}{2} (x_{0} - p) \times (x_{1} - p) \cdot n = {\begin{array}{r} \frac{1}{2} ‖ (x_{0} - p) \times (x_{1} - p) ‖ & inside \\ - \frac{1}{2} ‖ (x_{0} - p) \times (x_{1} - p) ‖ & outside \end{array} \end{matrix}

Signed Areas:

\begin{aligned} A_{2} = \frac{1}{2} (x_{0} - p) \times (x_{1} - p) \cdot n \\ A_{0} = \frac{1}{2} (x_{1} - p) \times (x_{2} - p) \cdot n \\ A_{1} = \frac{1}{2} (x_{2} - p) \times (x_{0} - p) \cdot n \\ A = A_{0} + A_{1} + A_{2} \end{aligned}

$\vb{p}$ $b_0 = A_0/A$ $b_1 = A_1/A$ $b_2 = A_2/A$ $b_0 + b_1 +b_2 = 1$ )

$\vb{p} = b_0\vb{x}_0 + b_1\vb{x}_1 + b_2 \vb{x}_2$

=> Gourand Shading: Using barycentric interpolation (no longer popular)

Example: Tetrahedral Volume

$\vb{x}_{10} = \vb{x}_1 - \vb{x}_0$ $\vb{x}_{20} = \vb{x}_2 - \vb{x}_0$ $\vb{x}_{20} = \vb{x}_2 - \vb{x}_0$
$A = \frac{1}{2} ||\vb{x}_{10} - \vb{x} _ {20}||$
$h = \vb{x}_{30}\cdot \vb{n} = \vb{x}_{30} \cdot \frac{\vb{x}_{10} - \vb{x} _ {20}}{||\vb{x}_{10} - \vb{x} _ {20}||}$
Volume: (signed)
$\begin{matrix} V = \frac{1}{3} h A = \frac{1}{6} x_{30} \cdot x_{10} \times x_{20} = \frac{1}{6} [\begin{matrix} x_{1} & x_{2} & x_{3} & x_{0} \\ 1 & 1 & 1 & 1 \end{matrix}] \end{matrix}$

Example: Barycentric Weight in Tetrahedra

$\vb{p}$ $V_0 = \mathrm{Vol}(\vb{x}_3, \vb{x}_2, \vb{x}_1, \vb{p})$ , ….

$\vb{p}$ $V_0, V_1,V_2,V_3 >0$

$b_n = V_n /V$ $\vb{p} = b_0\vb{x}_0 + b_1\vb{x}_1 + b_2\vb{x}_2 + b_3\vb{x}_3$

Example: Particle-Triangle Intersection

$t$ $(\vb{p}(t) - \vb{x}_0)\cdot \vb{x}_{10}\times\vb{x}_{20} = 0$ $\vb{p}(t) = \vb{p} + t\vb{v}$
$t = \frac{(p - x_{0}) \cdot x_{10} \times x_{20}}{v \cdot x_{10} \times x_{20}}$
$\vb{p}(t)$ inside or outside

Matrices

Definitions

\begin{matrix} A = [\begin{array}{lll} a_{00} & a_{01} & a_{02} \\ a_{10} & a_{11} & a_{12} \\ a_{20} & a_{21} & a_{22} \end{array}] = [\begin{array}{lll} a_{0} & a_{1} & a_{2} \end{array}] \in R^{3 \times 3} \end{matrix}

Transpose / Diagonal / Identity / Symmetric:
$\begin{matrix} A^{T} = [\begin{array}{lll} a_{00} & a_{01} & a_{02} \\ a_{10} & a_{11} & a_{12} \\ a_{20} & a_{21} & a_{22} \end{array}]; [\begin{matrix} a_{00} \\ a_{11} \\ a_{22} \end{matrix}]; I = [\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}]; A^{T} = A \end{matrix}$

Orthogonality

An orthogonal matrix is a matrix made of orthogonal unit vectors.

\begin{matrix} A = [\begin{matrix} a_{0} & a_{1} & a_{2} \end{matrix}], such that a_{i}^{T} a_{j} = {\begin{aligned} 1, & if i = j \\ 0, & if i \neq j \end{aligned} \end{matrix}

\begin{matrix} A^{T} A = [\begin{matrix} a_{0}^{T} \\ a_{1}^{T} \\ a_{2}^{T} \end{matrix}] [\begin{array}{lll} a_{0} & a_{1} & a_{2} \end{array}] = [\begin{array}{ccc} a_{0}^{T} a_{0} & a_{0}^{T} a_{1} & a_{0}^{T} a_{2} \\ a_{1}^{T} a_{0} & a_{1}^{T} a_{1} & a_{1}^{T} a_{2} \\ a_{2}^{T} a_{0} & a_{2}^{T} a_{1} & a_{2}^{T} a_{2} \end{array}] = I; A^{T} = A^{- 1} \end{matrix}

Matrix Transformation

Rotation can be represented by an orthogonal matrix
(can represent local -> world)
(Considering of local coord. vect.)
$\begin{matrix} \begin{aligned} A x = u \\ A y = v \\ A z = w \end{aligned}} \Rightarrow A = [\begin{matrix} u & v & w \end{matrix}] \end{matrix}$
Scaling can be represented by a diagonal matrix
(Consisting of scaling factors)
$\begin{matrix} D = [\begin{matrix} d_{x} \\ d_{y} \\ d_{z} \end{matrix}] \end{matrix}$
Singular Value Decomposition
$\vb{A} = \vb{UDV}^{\mathrm{T}}$ $\vb{D}$ $\vb{U}$ $\vb{V}$ - Orthogonal)
$\vb{V}^{\mathrm{T}}$ $\vb{D}$ ) -> Rotation (All can be decomposed as rotaion and scaling (even in 3D))
Eigenvalue Decomposition
symmetric $\vb{A} = \vb{UDU}^{-1} = \vb{UDU}^{\mathrm{T}}$
$\vb{U} = \begin{bmatrix}\cdots & \vb{u}_i & \cdots\end{bmatrix}$ $\vb{u}_i$ $d_i$
$\begin{matrix} A u_{i} = {U D U}^{T} u_{i} = U D [\begin{matrix} ⋮ \\ 0 \\ 1 \\ 0 \\ ⋮ \end{matrix}] = U [\begin{matrix} ⋮ \\ 0 \\ d_{i} \\ 0 \\ ⋮ \end{matrix}] = d_{i} u_{i} \end{matrix}$
For asymmetric matrices -> eigenvalue can be imaginary nums
Symmetric Positive Definiteness (s.p.d.)
-> Linear system
For a s.p.d., if only if all its eigenvalues are positive
- $\vb{A}$ $\vb{v}^{\mathrm{T}}\vb{Av} >0$ $\vb{v} \ne 0$
- $\vb{A}$ $\vb{v}^{\mathrm{T}}\vb{Av} \ge 0$ $\vb{v} \ne 0$
Meaning:
- $d>0$ $\vb{v}^{\mathrm{T}}d\vb{v} > 0$ $\vb{v} \ne0$ ;
- $d_0, d_1, ... >0$ $\vb{v} \ne0$ )
  $\begin{matrix} \Rightarrow v^{T} D v = v^{T} [\begin{matrix} ⋱ \\ d_{i} \\ ⋱ \end{matrix}] v = \sum d_{i} v_{i}^{2} > 0 \end{matrix}$
- $d_0, d_1, ... >0$ $\vb{U}$ $\vb{v} \ne0$ )
  $\Rightarrow v^{T} ({U D U}^{T}) v = v^{T} {U U}^{T} ({U D U}^{T}) {U U}^{T} v = (U^{T} v)^{T} (D) (U^{T} v) > 0$
diagonally domiant $a_{ii} > \sum_{i\ne j}|a_{ij}|$ $i$ )
- Properties:
- $\vb{A}$ $\vb{B} = \begin{bmatrix}\vb{A}&\vb{-A}\\ \vb{-A}& \vb{A}\end{bmatrix}$ is symmetric semi-definite.

Linear Solver

$\vb{Ax} = \vb{b}$ $\vb{A}$ $\vb{x}$ $\vb{b}$ - Boundary vector)

$\vb{A}^{-1}$ $\vb{A}$ $\vb{x} = \vb{A}^{-1}\vb{b}$ )

Direct Linear Solver

LU factorization (Alt: Choleskey, LDL^T, etc.)

\begin{matrix} A = L U = \underset{L}{\underset{⏟}{[\begin{array}{ccc} l_{00} \\ l_{10} & l_{11} \\ ⋮ & \dots & ⋱ \end{array}]}} \underset{U}{\underset{⏟}{[\begin{array}{ccc} ⋱ & \dots & ⋮ \\ u_{n - 1, n - 1} & u_{n - 1, n} \\ u_{n, n} \end{array}]}} \end{matrix}

First solve: (up -> down)

\begin{matrix} L y = b \Rightarrow [\begin{array}{ccc} l_{00} \\ l_{10} & l_{11} \\ ⋮ & \dots & ⋱ \end{array}] [\begin{matrix} y_{0} \\ y_{1} \\ ⋮ \end{matrix}] = [\begin{matrix} b_{0} \\ b_{1} \\ ⋮ \end{matrix}] \Rightarrow {\begin{aligned} y_{0} = b_{0} / l_{00} \\ y_{1} = (b_{1} - l_{10} y_{0}) / l_{11} \\ . . . \end{aligned} \end{matrix}

Then solve: (down -> up)

\begin{matrix} U x = y \Rightarrow [\begin{array}{ccc} ⋱ & \dots & ⋮ \\ u_{n - 1, n - 1} & u_{n - 1, n} \\ u_{n, n} \end{array}] [\begin{matrix} ⋮ \\ x_{n - 1} \\ x_{n} \end{matrix}] = [\begin{matrix} ⋮ \\ y_{n - 1} \\ y_{n} \end{matrix}] \Rightarrow {\begin{aligned} x_{n} = y_{n} / u_{n, n} \\ x_{n - 1} = (y_{n - 1} - u_{n - 1, n} x_{n}) / u_{n - 1, n - 1} \\ \dots \end{aligned} \end{matrix}

Properties:

$\vb{A}$ $\vb{L}$ $\vb{U}$ are not that sparse, dep on the permutation (modify the order) -> MATLAB (LUP)
$\vb{A}$ , factorization could be done once (Save costs)
Hard to parallelized (Intel MKL PARDISO)

Iterative Linear Solver

x^{[k + 1]} = x^{[k]} + α M^{- 1} (b - A x^{[k]})

$\alpha$ $\vb{M}$ $\vb{b - Ax}^{[k]}$ - Residual error (for perfect solution -> = 0))

$\vb{b-Ax}^{[0]} = \text{const}$ at first)

\begin{aligned} b - A x^{[k + 1]} & = b - A x^{[k]} - α A M^{- 1} (b - A x^{[k]}) \\ = (I - α A M^{- 1}) (b - A x^{[k]}) = {(I - α A M^{- 1})}^{k + 1} (b - A x^{[0]}) \end{aligned}

$\rho\left(\mathbf{I}-\alpha \mathbf{A M}^{-1}\right)$ - Spectral radius (largest absolute value of eigenvalues))

b - A x^{[k + 1]} \to 0, if ρ (I - α {A M}^{- 1}) < 1

$\vb{M}$ $\vb{M} = \mathrm{diag}(\vb{A})$ Jacobi $\vb{M} = \text{lower}(\vb{A})$ (Gauss-Seidel)

Accelerate converge methods: Chebyshev, Conjugate Gradient, …

Properties:

Simple; Fast for inexact sol; Parallelable
Convergence condition (not converge for every matrix); Slow for exact solutions

Tensor Calculus

Basic Concepts

1st-Order Derivatives
$f(\vb{x})\in \R$ , then
$\begin{matrix} d f = \frac{\partial f}{\partial x} d x + \frac{\partial f}{\partial y} d y + \frac{\partial f}{\partial z} d z = [\begin{matrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \end{matrix}] [\begin{array}{l} d x \\ d y \\ d z \end{array}]; \frac{\partial f}{\partial x} = [\begin{array}{lll} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \end{array}] or \nabla f (x) = [\begin{matrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z} \end{matrix}] \end{matrix}$
$\vb{f(x)} = \begin{bmatrix}f(\vb{x})\\ g(\vb{x}) \\ h(\vb{x})\end{bmatrix}\in \R^3$ , then
$\begin{matrix} Jacobian: J (x) = \frac{\partial f}{\partial x} = [\begin{array}{lll} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \end{array}]; Divergence: \nabla \cdot f = \frac{\partial f}{\partial x} + \frac{\partial g}{\partial y} + \frac{\partial h}{\partial z}; Curl: \nabla \times f = [\begin{array}{l} \frac{\partial h}{\partial y} - \frac{\partial g}{\partial z} \\ \frac{\partial f}{\partial z} - \frac{\partial h}{\partial x} \\ \frac{\partial g}{\partial x} - \frac{\partial f}{\partial y} \end{array}] \end{matrix}$
2nd-Order Derivatives
$f(\vb{x})\in \R$ , then (Hessian is symmetric, tangent stiffness)
$\begin{matrix} Hessian: H = J (\nabla f (x)) = [\begin{array}{ccc} \frac{\partial^{2} f}{\partial x^{2}} & \frac{\partial^{2} f}{\partial x \partial y} & \frac{\partial^{2} f}{\partial x \partial z} \\ \frac{\partial^{2} f}{\partial x \partial y} & \frac{\partial^{2} f}{\partial y^{2}} & \frac{\partial^{2} f}{\partial y \partial z} \\ \frac{\partial^{2} f}{\partial x \partial z} & \frac{\partial^{2} f}{\partial y \partial z} & \frac{\partial^{2} f}{\partial z^{2}} \end{array}]; Laplacian: \nabla \cdot \nabla f (x) = \nabla^{2} f (x) = Δ f (x) = \frac{\partial^{2} f}{\partial x^{2}} + \frac{\partial^{2} f}{\partial y^{2}} + \frac{\partial^{2} f}{\partial z^{2}} \end{matrix}$
Taylor Expansion
$f({x})\in \R$ , then
$f (x) = f (x_{0}) + \frac{\partial f (x_{0})}{\partial x} (x - x_{0}) + \frac{1}{2} \frac{\partial f^{2} (x_{0})}{\partial x^{2}} {(x - x_{0})}^{2} + \dots$
$f(\vb{x})\in \R$ , then
$\begin{aligned} f (x) & = f (x_{0}) + \frac{\partial f (x_{0})}{\partial x} (x - x_{0}) + \frac{1}{2} {(x - x_{0})}^{T} \frac{\partial f^{2} (x_{0})}{\partial x^{2}} (x - x_{0}) + \dots \\ = f (x_{0}) + \nabla f (x_{0}) \cdot (x - x_{0}) + \frac{1}{2} {(x - x_{0})}^{T} H (x - x_{0}) + \dots \end{aligned}$
$\vb{H}$ is s.p.d., second order derivative > 0 => interesting properties (to be discussed)

Lecture 3 Rigid Body Dynamics

(Single rigid body: dynamics / rotation / …)

Rigid bodies: Assume no deformations

$\vb{s}^{[k]}$ ove

Translation Motion

$\vb{x}$ $\vb{v}$ $M$ - Mass; Force can be the function of pos, vel, t, …) -> Solve integral

\begin{matrix} {\begin{aligned} v (t^{[1]}) = v (t^{[0]}) + M^{- 1} \int_{t^{[0]}}^{t^{[1]}} f (x (t), v (t), t) d t \\ x (t^{[1]}) = x (t^{[0]}) + \int_{t^{[0]}}^{t^{[1]}} v (t) d t \end{aligned} \end{matrix}

Integration Methods Explained

$\vb{x}(t) = \int \vb{v}(t)\ \mathrm d t$ is the area.

Explicit Euler $t^{[0]}$
$\int_{t^{[0]}}^{t^{[1]}} v (t) d t \approx Δ t v (t^{[0]})$
$\int_{t^{[0]}}^{t^{[1]}} \mathbf{v}(t) d t=\Delta t \mathbf{v}\left(t^{[0]}\right)+\frac{\Delta t^{2}}{2} \mathbf{v}^{\prime}\left(t^{[0]}\right)+\cdots = \Delta t \mathbf{v}\left(t^{[0]}\right)+ \mathcal O(\Delta t^2)$ $\mathcal O(\Delta t^2)$ )
Implicit Euler $t^{[1]}$
$\int_{t^{[0]}}^{t^{[1]}} v (t) d t \approx Δ t v (t^{[1]})$
$\int_{t^{[0]}}^{t^{[1]}} \mathbf{v}(t) d t=\Delta t \mathbf{v}\left(t^{[1]}\right)+\frac{\Delta t^{2}}{2} \mathbf{v}^{\prime}\left(t^{[1]}\right)+\cdots = \Delta t \mathbf{v}\left(t^{[1]}\right)+ \mathcal O(\Delta t^2)$ $\mathcal O(\Delta t^2)$ )
Mid-point $t^{[0.5]}$
$\int_{t^{[0]}}^{t^{[1]}} v (t) d t \approx Δ t v (t^{[0.5]})$
Taylor:
$\begin{aligned} \int_{t^{[0]}}^{t^{[1]}} v (t) d t & = \int_{t^{[0]}}^{t^{[0.5]}} v (t) d t + \int_{t^{[0.5]}}^{t^{[1]}} v (t) d t \\ = \frac{1}{2} Δ t v (t^{[0.5]}) - \frac{Δ t^{2}}{2} v^{' (t^{[0.5]})} + O (Δ t^{3}) + \frac{1}{2} Δ t v (t^{[0.5]}) + \frac{Δ t^{2}}{2} v^{' (t^{[0.5]})} + O (Δ t^{3}) \\ = Δ t v (t^{[0.5]}) + O (Δ t^{3}) \end{aligned}$
Final Method in this case: Semi-implicit (Mid-point)
Velocity -> Explicit; Position -> Implicit
$\begin{matrix} {\begin{aligned} v^{[1]} & = v^{[0]} + Δ t M^{- 1} f^{[0]} \\ x^{[1]} & = x^{[0]} + Δ t v^{[1]} \end{aligned} \end{matrix}$
Alternative: Leapfrog Integration
v and x not overlap (Mid-points)
$\begin{matrix} {\begin{aligned} v^{[0.5]} & = v^{[- 0.5]} + Δ t M^{- 1} f^{[0]} \\ x^{[1]} & = x^{[0]} + Δ t v^{[0.5]} \end{aligned} \end{matrix}$

Type of Forces

Gravity Force $\vb{f}_{\text{gravity}}^{[0]} = M\vb{g}$
$\vb{f}^{[0]}_{\text{drag}} = -\sigma \vb{v}^{[0]}$ $\sigma$ - drag coefficient) -> Reduced by following
$\vb{v}^{[1]} = \alpha \vb{v}^{[0]}$ $\alpha$ - decay coefficient)

Translation Only Simulation

Steps $M$ $\Delta t$ are user spec var)

$\vb{f}_i^{[0]} \leftarrow \text{Force} (\vb{x}_i^{[0]}, \vb{v}_i^{[0]})$
$\vb{f}^{[0]} \leftarrow \sum \vb{f}_i^{[0]}$
$\vb{v} ^{[1]} \leftarrow \vb{v}^{[0]} + \Delta tM^{-1}\vb{f}^{[0]}$
$\vb{x}^{[1]} \leftarrow \vb{x}^{[0]} + \Delta t\vb{v}^{[1]}$

Rotation Motion

Rotation Representations

Rotation Represented by Matrix

\begin{matrix} R = [\begin{array}{lll} r_{00} & r_{01} & r_{02} \\ r_{10} & r_{11} & r_{12} \\ r_{20} & r_{21} & r_{22} \end{array}] \end{matrix}

Suitable in graphics, rotation for vertices

Not suitable for dynamics:

Too much redundancy: 9 elem, 3 dof
Not intuitive
Defining its time derivative (rotational vel) is difficult

Rotation Represented by Euler Angles

Use 3 axial rotations to represent one general rotation. Each axial rotation uses an angle.

Used in Unity. (the order is rotaion-by-Z / X / Y) Intuitive.

Not suitable for dynamics:

Lose DOFs in certain statues: Gimbal lock
Defining its time derivative is difficult

Rotation Represented by Quaternion

Complex multiplications: In the complex system, two numbers represent a 2D point. => Quaternion: i, j, k are imaginary numbers (3D space) Four numbers represent a 3D point (with multiplication and division).

Arithematic

$\vb{q} = [s\ \ \vb{v}]$ $s$ $\vb{v}$ $\vb{ijk}$

$a \mathbf{q}=\left[\begin{array}{ll} a s & a \mathbf{v} \end{array}\right]$ Scalar-quaternion Multiplication
$\mathbf{q}_{1} \pm \mathbf{q}_{2}=\left[s_{1} \pm s_{2} \quad \mathbf{v}_{1} \pm \mathbf{v}_{2}\right]$ Addition/Subtraction (Same as vector)
$\mathbf{q}_{1} \times \mathbf{q}_{2}=\left[s_{1} s_{2}-\mathbf{v}_{1} \cdot \mathbf{v}_{2} \quad s_{1} \mathbf{v}_{2}+s_{2} \mathbf{v}_{1}+\mathbf{v}_{1} \times \mathbf{v}_{2}\right]$ Multiplication
$||\vb{q}|| = \sqrt{s^2 + \vb{v}\cdot\vb{v}}$ Magnitude

In Unity: provide multiplication, but no addition/subtraction/…; Use w, x, y, z -> s, v

Representation

$\vb{v}$ $\theta$

\begin{matrix} {\begin{cases} q = [\begin{array}{ll} \cos \frac{θ}{2} & v \end{array}] \\ ∥ q ∥ = 1 \end{cases} \Rightarrow {\begin{cases} q = [\begin{array}{ll} \cos \frac{θ}{2} & v \end{array}] \\ ∥ v ∥^{2} = \sin^{2} \frac{θ}{2} \end{cases} \end{matrix}

Convertible to the matrix:

\begin{matrix} R = [\begin{array}{ccc} s^{2} + x^{2} - y^{2} - z^{2} & 2 (x y - s z) & 2 (x z + s y) \\ 2 (x y + s z) & s^{2} - x^{2} + y^{2} - z^{2} & 2 (y z - s x) \\ 2 (x z - s y) & 2 (y z + s x) & s^{2} - x^{2} - y^{2} + z^{2} \end{array}] \end{matrix}

Rotation Representations in Unity

Rotation Motion

$\mathrm{\boldsymbol {\omega}}$ to denote angular velocity:

$\boldsymbol {\omega}$ $\boldsymbol {\omega}$ -> the speed (sim to the representation of quaternion)

Torque and Inertia

Torque

$\vb{r}_i$ $\vb{Rr}_i$ $\vb{f}_i$ is a force)

The rotational equiv of a force, describing the rotational tendency caused by a force.

$\boldsymbol \tau_i$ $\vb{Rr}_i$ $\vb{f}_i$ $||\vb{Rr}_i||$ $||\vb{f}_i||$ $\sin \theta$ (the angle of the two vectors)

τ_{i} \leftarrow R r_{i} \times f_{i}

Inertia

Describes the resistance to rotational tendency caused by torque (not const)

Left side (heavier point far away from the torque) has higher resistance (inertia) to the rotational tendency, slower rotation

Ref state inertia, change with rotation (dep on pose). But no need to re-compute every time

I_{r e f} = \sum m_{i} (r_{i}^{T} r_{i} 1 - r_{i} r_{i}^{T})

$\vb{1}$ - 3x3 identity matrix)

\begin{aligned} I & = \sum m_{i} (r_{i}^{T} R^{T} R r_{i} 1 - R r_{i} r_{i}^{T} R^{T}) \\ = \sum m_{i} (R r_{i}^{T} r_{i} 1 R^{T} - R r_{i} r_{i}^{T} R^{T}) \\ = \sum m_{i} R (r_{i}^{T} r_{i} 1 - r_{i} r_{i}^{T}) R^{T} \\ = R I_{r e f} R^{T} \end{aligned}

Use torque to represent

Torque $\boldsymbol{\tau}_i = (\vb{Rr}_i)\times \vb{f}_i$
$\boldsymbol{\tau}_{i} = \sum\boldsymbol{\tau}_{i}$
The rotational equivalent of massinertia $\vb{I}$ (Result: 3x3):
$\vb{I}_{\mathrm{ref}} = \sum m_i(\vb{r}_i^{\mathrm{T}}\vb{r}_i\vb{1} - \vb{r}_i\vb{r}_i^{\mathrm{T}})$
$\vb{I} = \vb{R} \vb{I}_{\mathrm{ref}} \vb{R}^{\mathrm{T}}$

Rigid Body Simulation

Translational and Rotational Motion

Translation (Linear)
$\vb{v}$ $\vb{x}$ (transform.position in Unity)
$M$ $\vb{f}$
$\begin{matrix} {\begin{cases} v^{[1]} = v^{[0]} + Δ t M^{- 1} f^{[0]} \\ x^{[1]} = x^{[0]} + Δ t v^{[1]} \end{cases} \end{matrix}$
$||\vb{q}||\ne1$ (in Unity automatically)
$\boldsymbol{\omega}$ $\vb{q}$ (transform.rotation in Unity)
$\vb{I}$ $\boldsymbol{\tau}$
$\begin{matrix} {\begin{cases} ω^{[1]} = ω^{[0]} + Δ t {(I^{[0]})}^{- 1} τ^{[0]} \\ q^{[1]} = q^{[0]} + [0 \frac{Δ t}{2} ω^{[1]}] \times q^{[0]} \end{cases} \end{matrix}$

Rigid Body Simulation Process

In Unity: No 3x3 matrices, only 4x4 (use 4x4 and set the last col / line ); Provide .inverse to inverse; …

Implementation

$\vb{s = \{ v,x},\boldsymbol{\omega},\vb{q\}}$

Issues

Translation is easier, code translation first
$\boldsymbol{\omega}$ $\vb{q}$ , in this case the object will spin constantly
Gravity does NOT cause torque (except for air drag force)

Lecture 4 Rigid Body Contacts

Particle Collision Detection and Response

Distance Functions

Signed Distance Function

$\phi(\vb{x})$ to define the distance indicating which side as well (corresponding to 0 surface)

Examples

ϕ (x) = (x - p) \cdot n ϕ (x) = | | x - c | | - r ϕ (x) = \sqrt{| | x - p | |^{2} - ((x - p) \cdot n)^{2}} - r

Intersection of Signed Distance Functions

=> Bool operations

$\phi_0(\vb{x})< 0$ $\phi_1(\vb{x})< 0$ $\phi_2(\vb{x})< 0$ $\phi(\vb{x})=\max(\phi_0{\vb{x}}, \phi_1(\vb{x}), \phi_2(\vb{x}))$ $\phi(\vb{x}) = ?$ // not relavent (no collision)

Union of Signed Distance Functions

$\phi_0(\vb{x})< 0$ $\phi_1(\vb{x})< 0$ $\phi(\vb{x}) \approx \min(\phi_0(\vb{x}),\phi_1(\vb{x}))$ // approximate -> correct near outer boundary $\phi(\vb{x}) = \min(\phi_0(\vb{x}), \phi_1(\vb{x}))$

-> We can consider collision detection with the union of two objects as collision detection with two separate objects

Penalty methods

(Implicit integration is better)

Quadratic Penalty Method

Check if collide - Yes -> Apply a force at the point (in the next update)

For quadratic penalty (strength) potential, the force is linear

Problem: Already inside -> cause artifacts

=> Add buffer help less the penetration issue (cannot strictly prevent)

$k$ $k$ too high -> too much force generated (overshooting)

Log-Barrier Penalty Method

$\phi(\vb{x})< 0$ $\Delta t$

$\vb{f} \leftarrow \rho \frac{1}{\phi(\vb{x})} \vb{N}$ $\rho$ - Barrier strength)

Problems: cannot prevent overshooting when very close to the surface; when penatration occurs, the penatration will be higher and higher (-> smaller step size -> higher costs)

=> Log-Barrier limited within a buffer as well to solve

Frictional contacts are difficult to handle

Impulse method

Update the vel and pos as the collision occurs

Changing the position:

velocity $\mu_{\vb{N}}$ $\in [0,1]$ $a$ - frictional decay of vel)

$a$ should be minimized but not violating Coulomb’s law

\begin{matrix} ‖ v_{T}^{new} - v_{T} ‖ \leq μ_{T} ‖ v_{N}^{new} - v_{N} ‖ \\ (1 - a) ‖ v_{T} ‖ \leq μ_{T} (1 + μ_{N}) ‖ v_{N} ‖ \end{matrix}

Therefore:

a ⟵ max (1 - \underset{dynamic friction}{\underset{⏟}{μ_{T} (1 + μ_{N}) ‖ v_{N} ‖ / ‖ v_{T} ‖}}, \underset{static friction}{\underset{⏟}{0}})

Can precisely control the friction effects

Rigid Collision Detection and Response by Impulse

Rigid Body Collision Detection

$\vb{x}_i \leftarrow \vb{x}+ \vb{Rr}_i$ (from the mass center to the vertices)

$\phi(\vb{x})<0$

Rigid Body Collision Response by Impulse

$i$ $\vb{v}$ $\boldsymbol{\omega}$ - angular vel)

\begin{matrix} {\begin{aligned} x_{i} & \leftarrow x + R r_{i} & (Position) \\ v_{i} & \leftarrow v + ω \times R r_{i} & (Velocity) \end{aligned} \end{matrix}

$\vb{x}_i$ $\vb{v}_i$ directly since they are not state var. (in Lec.3: 4 var are vel of mass center / pos of mass center / rotational pose / angular vel)

Applying an impulse $\vb{j}$ $i$ $\Delta \vb{v} = \frac{\vb{j}}{M}$ $\vb{Rr}_i \times \vb{j}$ $\vb{j}$ )

\begin{matrix} {\begin{aligned} v^{new} & \leftarrow v + \frac{1}{M} j \\ ω^{new} & \leftarrow ω + I^{- 1} (R r_{i} \times j) \end{aligned} \end{matrix}

\begin{aligned} \Rightarrow v_{i}^{new} & = v^{new} + ω^{new} \times R r_{i} \\ = v + \frac{1}{M} j + (ω + I^{- 1} (R r_{i} \times j)) \times R r_{i} \\ = v_{i} + \frac{1}{M} j - (R r_{i}) \times (I^{- 1} (R r_{i} \times j)) \end{aligned}

Cross Product as a Matrix Product

$\vb{r}\times$ $\vb{r}^*$

\begin{matrix} r \times q = [\begin{array}{l} r_{y} q_{z} - r_{z} q_{y} \\ r_{z} q_{x} - r_{x} q_{z} \\ r_{x} q_{y} - r_{y} q_{x} \end{array}] = [\begin{array}{ccc} 0 & - r_{z} & r_{y} \\ r_{z} & 0 & - r_{x} \\ - r_{y} & r_{x} & 0 \end{array}] [\begin{array}{l} q_{x} \\ q_{y} \\ q_{z} \end{array}] = r^{*} q \end{matrix}

In our case:

\Rightarrow v_{i}^{new} = v_{i}^{new} = v_{i} + \frac{1}{M} j - {(R r_{i})}^{*} I^{- 1} {(R r_{i})}^{*} j

$\vb{K}$ $\vb{j}$ can be computed with the following equations.

\begin{matrix} \Rightarrow v_{i}^{new} - v_{i} = K j \\ K \leftarrow \frac{1}{M} 1 - {(R r_{i})}^{*} I^{- 1} {(R r_{i})}^{*} \end{matrix}

Implementation

Other details:

For several vertices in collision, use their average
decrease the restitution $\mu_{\vb{N}}$ to reduce oscillation
Don’t update the position: not linear

Shape Matching

Basic Idea

Allow each vertex to have its own velocity, so it can move by itself

Move vertices independently by its velocity, with collision and friction being handled (use the impulse method for every point)
Enforce the rigidity constraint to become a rigid body again (IMPORTANT)

Mathematical Formulation

$\vb{c}$ $\vb{R}$ $\vb{c}$ - mass center): want the final rigid body (a square) close enough to the trapozoid

$\vb{A}$ $\sum \vb{Ar}_i = 0$ $E$ $= \frac{1}{2}\left\|\mathbf{c}+\mathbf{A} \mathbf{r}_{i}-\mathbf{y}_{i}\right\|^{2}$ )

\begin{aligned} {c, R} & = argmin \sum_{i} \frac{1}{2} {‖ c + R r_{i} - y_{i} ‖}^{2} \\ \Rightarrow {c, A} & = argmin \sum_{i} \frac{1}{2} {‖ c + A r_{i} - y_{i} ‖}^{2} \end{aligned}

$\vb{c}$ $\vb{A}$ (Find derivatives):

\begin{matrix} \frac{\partial E}{\partial c} = \sum_{i} c + A r_{i} - y_{i} = \sum_{i} c - y_{i} = 0 \\ \Rightarrow c = \frac{1}{N} \sum_{i} y_{i} (average) \end{matrix}

\begin{matrix} \frac{\partial E}{\partial A} = \sum_{i} (c + A r_{i} - y_{i}) r_{i}^{T} = 0 \\ A = (\sum_{i} (y_{i} - c) r_{i}^{T}) {(\sum_{i} r_{i} r_{i}^{T})}^{- 1} \overset{Polar Decomposition}{=} \underset{rotation}{\underset{⏟}{R}} \underset{deformation}{\underset{⏟}{S}} \end{matrix}

Remind $\vb{A} = \vb{UDV}^{\mathrm{T}}$ (rotaion, scaling and rotation)

$\vb{A} = (\vb{UV}^{\mathrm{T}})(\vb{VDV}^{\mathrm{T}}) = \vb{RS}$ $\vb{VDV}^{\mathrm{T}} = \vb{S}$ )

Implementation

Physical quantities are attached to each vertex, not to the entire body.

$\mathrm{Polar}(\vb A)$ is provided, ultilizing the polar deposition technique)

Properties:

Easy to implement and compatible with other nodal systems: cloth, soft bodies, particle fluids, …
Difficult to strictly enforce friction and other goals. The rigidification process will destroy them (may require iter)
More suitable when the friction accuracy is unimportant, i.e., buttons on clothes

Lecture 5 Physics-Based Cloth Simulation

A Mass-Spring System

Spring Systems

An Ideal Spring

$k$ - Spring Stiffness)

E (x) = \frac{1}{2} k (x - L)^{2}; f (x) = - \frac{d E}{d x} = - k (x - L)

\begin{matrix} E (x) = \frac{1}{2} k (∥ x_{i} - x_{j} ∥ - L) \\ f_{i} = - f_{j} {\begin{array}{r} f_{i} (x) = - \nabla_{i} E = - k (‖ x_{i} - x_{j} ‖ - L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖} \\ f_{j} (x) = - \nabla_{j} E = - k (‖ x_{j} - x_{i} ‖ - L) \frac{x_{j} - x_{i}}{‖ x_{j} - i_{j} ‖} \end{array} \end{matrix}

Multiple Springs

The energies and forces can be simply summed up

\begin{matrix} E = \sum_{j = 0}^{3} E_{j} = \sum_{j = 0}^{3} (\frac{1}{2} k {(‖ x_{i} - x_{j} ‖ - L)}^{2}) \\ f_{i} = - \nabla_{i} E = \sum_{j = 0}^{3} (- k (‖ x_{i} - x_{j} ‖ - L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖}) \end{matrix}

Structures in Simulations

Structured Spring Networks

Unstructured Spring Networks

Unstructured triangle mesh

-> the edges into spring networks (usually in cloth simulations)

Blue lines for bending resistance (every neighboring triangle pair)

Triangle Mesh Representation

Two arrays: Vertex & Triangle lines

$\vb{x}_0$ $\vb{x}_1$ $\vb{x}_2$ $\vb{x}_3$ $\vb{x}_4$ } (3D vectors)
Triangle list: {1, 2, 3, 0, 1, 3, 0, 3, 4} (Index triples) => ({1, 2, 3} for Triangle 0; {0, 1, 3} for Triangle 1; …)

Not only edges but also an inner one (each triangle has 3 edges but there are repeated ones)

Topological Construction

Sort triangle edge triples: edge vertex index 0 & 1 and triangle index (index 0 < index)

Triple list:


xxxxxxxxxx
{{1, 2, 0}, {2, 3, 0}, {1, 3, 0},   // Green triangle: {1, 2} for edge index; {0} at the end for triangle index
 {0, 1, 1}, {1, 3, 1}, {0, 3, 1},
 {0, 3, 2}, {3, 4, 2}, {0, 4, 2}}

Sorted triple list: {{0, 1, 1}, {0, 3, 1}, {0, 3, 2}, {0, 4, 2}, {1, 2, 0}, {1, 3, 0}, {1, 3, 1}, {2, 3, 0}, {3, 4, 2}}
Repeated edges {1, 3, 0}|{1, 3, 1} & {0, 3, 1}|{0, 3, 2} shows in the neighbor => eliminate
Final edge list: {{0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 3}, {2, 3}, {3, 4}}
Neighboring triangle list: {{1, 2}, {0, 1}} (for bending) (or use neighboring edge list)

Integrators

Explicit Integration

Scheme

$m_i$ $i$ $E$ $L$ - Edge length list (pre-computed))

For every vertex:
- $\vb{f}_i \leftarrow \text{Force}(\vb{x}_i, \vb{v}_i)$
- $\vb{v}_i \leftarrow \vb{v}_i + \Delta tm_i^{-1}\vb{f}_i$
- $\vb{x}_i \leftarrow \vb{x}_i + \Delta t\vb{v}_i$
$e$ )
- $i\leftarrow E[e][0]$ $[e]$ $[0]$ )
- $j \leftarrow E[e][1]$
- $L_e \leftarrow L[e]$
- $\vb{f}\leftarrow -k(\|\vb{x}_i - \vb{x}_j\| - L_e) \frac{\vb{x}_i - \vb{x}_j}{\|\vb{x}_i - \vb{x}_j\|}$
- $\vb{f}_i\leftarrow \vb{f}_i + \vb{f}$
- $\vb{f}_j\leftarrow \vb{f}_j + \vb{f}$

Problem

$k$ $\Delta t$ is too large

$\Delta t$ => Slow down the simulation

Implicit Integration

General Scheme (Euler Method)

$\vb{x}$ $\vb{v}$ $\vb{M}\in \R^{3N \times 3N}$ - Mass matrix (usually diagonal))

\begin{matrix} {\begin{cases} v^{[1]} = v^{[0]} + Δ t M^{- 1} f^{[1]} \\ x^{[1]} = x^{[0]} + Δ t v^{[1]} \end{cases} or {\begin{cases} x^{[1]} = x^{[0]} + Δ t v^{[0]} + Δ t^{2} M^{- 1} f^{[1]} \\ v^{[1]} = (x^{[1]} - x^{[0]}) / Δ t \end{cases} \end{matrix}

$\vb{v}^{[1]}$ $\vb{x}^{[1]}$ are unknown for the current time step => Find the x and v:

$\vb{f}$ $\vb{x}$ $\vb{f}$ may not be a linear function)

x^{[1]} = x^{[0]} + Δ t v^{[0]} + Δ t^{2} M^{- 1} f (x^{[1]})

$\|\vb{x}\|^2_\vb{M} = \vb{x}^{\mathrm{T}}\vb{Mx}$ ) => Optimization prob => Numerical schemes (Usually only conservative forces can use this energy)

x^{[1]} = argmin F (x) for F (x) = \frac{1}{2 Δ t^{2}} {‖ x - x^{[0]} - Δ t v^{[0]} ‖}_{M}^{2} + E (x)

$F(\vb{x})$ reaches 0 for the min pt.)

\nabla F (x^{[1]}) = \frac{1}{Δ t^{2}} M (x^{[1]} - x^{[0]} - Δ t v^{[0]}) - f (x^{[1]}) = 0 \Rightarrow x^{[1]} - x^{[0]} - Δ t v^{[0]} - Δ t^{2} M^{- 1} f (x^{[1]}) = 0

The Optimization Problem

Newton-Raphson Method

${x}^{[1]} = \operatorname{argmin} F(x)$ $F(x)$ is Lipschitz continuous)

$x^{(k)}$ $0 = F'(x) \approx F'(x^{(k)}) + F''(x^{(k)})(x-x^{(k)})$ (Taylor Expansion)

$\vb{0} = \grad F(\vb{x}) \approx \grad F(\vb{x}^{(k)}) + \frac{\partial F^2 (\vb{x}^{(k)})}{\partial \vb{x}^2} (\vb{x}-\vb{x}^{(k)})$ )

Steps:

$x^{(0)}$
$k = 0\ ...\ k$
$\Delta x \leftarrow -(F''(x^{(k)}))^{-1} F'(x^{(k)})$ $\Delta \mathbf{x} \leftarrow-\left(\frac{\partial F^{2}\left(\mathbf{x}^{(k)}\right)}{\partial \mathbf{x}^{2}}\right)^{-1} \grad F\left(\mathbf{x}^{(k)}\right)$ )
$x^{(k+1)} \leftarrow x^{(k)} + \Delta x$
$|\Delta x|$ is small then break
$\vb{x}^{[1]} \leftarrow \vb{x}^{(k+1)}$

$F''(x^*) > 0$ $F''(x^*)<0$ )

$F''(x)>0$ $F(x)$ has only one minimum

Simulation by Newton’s Method

$F(\mathbf{x})=\frac{1}{2 \Delta t^{2}}\left\|\mathbf{x}-\mathbf{x}^{[0]}-\Delta t \mathbf{v}^{[0]}\right\|_{M}^{2}+E(\mathbf{x})$

Derivative:
$\nabla F (x^{(k)}) = \frac{1}{Δ t^{2}} M (x^{(k)} - x^{[0]} - Δ t v^{[0]}) - f (x^{(k)})$
$\vb{H}(\vb{x}^{(k)})$ - Energy Hessian)
$\frac{\partial^{2} F (x^{(k)})}{\partial x^{2}} = \frac{1}{Δ t^{2}} M + H (x^{(k)})$

Steps:

$\vb{x}^{(0)}$ $\vb{x}^{[0]}$ $\vb{x}^{[0] }+\Delta t\vb{v}^{[0]}$
$k = 0 \ ...\ k$ :
$\left(\frac{1}{\Delta t^{2}} \mathbf{M}+\mathbf{H}\left(\mathbf{x}^{(k)}\right)\right) \mathbf{x}=\left(-\frac{1}{\Delta t^{2}} \mathbf{M}\left(\mathbf{x}^{(k)}-\mathbf{x}^{[0]}- \Delta t \mathbf{v}^{[0]}\right)+\mathbf{f}\left(\mathbf{x}^{(k)}\right)\right.$
$\vb{x}^{(k+1)} \leftarrow \vb{x}^{(k)} + \Delta \vb{x}$
$\|\Delta \vb{x}\|$ is small then break
$\vb{x}^{[1]}\leftarrow \vb{x}^{(k+1)}$
$\vb{v}^{[1]}\leftarrow (\vb{x}^{[1]} - \vb{x}^{[0]})/ \Delta t$

Spring Hessian

$\partial^2 F / \partial \vb{x}^2$ ) => if p.d. so that the ONLY min point is found and has no maximum (sufficient but not neccesary condition)

\begin{matrix} H (x) = \sum_{e = {i, j}} [\begin{array}{cc} \frac{\partial^{2} E}{\partial x_{i}^{2}} & \frac{\partial^{2} E}{\partial x_{i} \partial x_{j}} \\ \frac{\partial^{2} E}{\partial x_{i} \partial x_{j}} & \frac{\partial^{2} E}{\partial x_{j}^{2}} \end{array}] = \sum_{e = {i, j}} [\begin{array}{cc} H_{e} & - H_{e} \\ - H_{e} & H_{e} \end{array}] \end{matrix}

$\vb{H}_e$ $(i,i)$ $-\vb{H}_e$ $(i,j)$ , …

$\vb{H}_e$ (3 x 3): @ Lec.2, P48

H_{e} = k \frac{x_{i j} x_{i j}^{T}}{{‖ x_{i j} ‖}^{2}} + k (1 - \frac{L}{‖ x_{i j} ‖}) (I - \frac{x_{i j} x_{i j}^{T}}{{‖ x_{i j} ‖}^{2}}); x_{i j} = x_{i} - x_{j}

$k \frac{\mathbf{x}_{i j} \mathbf{x}_{i j}^{\mathrm{T}}}{\left\|\mathbf{x}_{i j}\right\|^{2}}$ $\left(\mathbf{I}-\frac{\mathbf{x}_{i j} \mathbf{x}_{i j}^{\mathrm{T}}}{\left\|\mathbf{x}_{i j}\right\|^{2}}\right)$ $\vb{v}^{\mathrm{T}}$ $\vb{v}$ at both sides)

$\left(1-\frac{L}{\left\|\mathbf{x}_{i j}\right\|}\right)$ $\|\vb{x}_{ij}\| < L_e$ (when compressed)

Conclusion $\vb{H}_e$ $\vb{H}_e$ $\vb {A}$ may not be s.p.d. either

\begin{matrix} A = \frac{1}{Δ t^{2}} M + H (x) = \underset{s.p.d.}{\underset{⏟}{\frac{1}{Δ t^{2}} M}} + \sum_{e = {i, j}} \underset{may not be s.p.d.}{\underset{⏟}{[\begin{array}{ccc} ⋱ & ⋮ & ⋮ \\ H_{e} & - H_{e} \\ - H_{e} & H_{e} \\ ⋱ \end{array}]}} \end{matrix}

$\Delta t$ => more p.d., actually sim to explicit integration with smaller time step => more stable)

Positive Definiteness of Hessian

When a spring is compressed, the spring Hessian may not be positive definite => Multiple minima

$E(x) = \frac{1}{2}k(x-L)^2$ $E''(x) = k>0$

Enforcement of P.D.

$\vb{A}$ $\vb{A}\Delta \vb{x} = \vb{b}$

Solution:

$\|\vb{x}_{ij} \|<L_e$ :
$H_{e} = k \frac{x_{i j} x_{i j}^{T}}{{‖ x_{i j} ‖}^{2}} + k (1 - \frac{L}{‖ x_{i j} ‖}) (I - \frac{x_{i j} x_{i j}^{T}}{{‖ x_{i j} ‖}^{2}})$
Choi and Ko. 2002. Stable But Responive Cloth. TOG (SIGGRAPH)

Linear Solvers

$\vb{A}\Delta \vb{x} = \vb{b}$

Jocobi Method

$\Delta \vb{x} \leftarrow 0$
$k = 0 \ ... \ k$ :
$\vb{r} \leftarrow \vb{b} - \vb{A}\Delta \vb{x}$ // Residual error
$\|r\|< \varepsilon$ $\varepsilon$
$\Delta \vb{x }\leftarrow \Delta \vb{x }+ \alpha \vb{D}^{-1}\vb{r}$ $\vb{D}$ $\vb{A}$

$\alpha = 1$ $\vb{A}$ : Diagonal Dominant

Other Solvers

Direct solvers (LU / LDLT / Cholesky / …) – [Intel MKL PARDISO]
- One shot, expensive but worthy if need exact sol
- $\vb{A}$
- Mostly suitable on CPUs
Iterative solvers
- Expensive to solve exactly, but controllable
- $\vb{A}$ , typically positive definiteness
- Suitable on both CPUs and GPUs
- Easy to implement
- Accelerable: Chebyshev, Nesterov, Conjugate Gradient …

After-Class Reading

Baraff and Witkin. 1998. Large Step in Cloth Simulation. SIGGRAPH.
One of the first papers using implicit integration. The paper proposes to use only one Newton iteration, i.e., solving only one linear system. This practice is fast, but can fail to converge.

Bending and Locking Issues

The Bending Spring Issue

A bending spring offers little resistance when cloth is nearly planar, since its length barely changes.

A Dihedral Angle Model

$\theta$ $\vb{f}_i = f(\theta) \vb{u}_i$

$\vb{u}_1$ $\vb{u}_2$ $\vb{n}_1$ $\vb{n}_2$
$\vb{u}_4 - \vb{u}_3$ $\vb{n}_1$ $\vb{n}_2$ )
$\vb{u}_1 + \vb{u}_2 + \vb{u}_3 + \vb{u}_4 = 0$ $\vb{u}_3$ $\vb{u}_4$ $\vb{n}_1$ $\vb{n}_2$

Conclusion:

$\vb{u}_1 = \|\vb{E}\|\frac{\vb{N}_1}{\|\vb{N}_1\|^2}$
$\vb{u}_2 = \|\vb{E}\|\frac{\vb{N}_2}{\|\vb{N}_2\|^2}$
$\mathbf{u}_{3}=\frac{\left(\mathbf{x}_{1}-\mathbf{x}_{4}\right) \cdot \mathbf{E}}{\|\mathbf{E}\|} \frac{\mathbf{N}_{1}}{\left\|\mathbf{N}_{1}\right\|^{2}}+\frac{\left(\mathbf{x}_{2}-\mathbf{x}_{4}\right) \cdot \mathbf{E}}{\|\mathbf{E}\|} \frac{\mathbf{N}_{2}}{\left\|\mathbf{N}_{2}\right\|^{2}}$
$\mathbf{u}_{4}=-\frac{\left(\mathbf{x}_{1}-\mathbf{x}_{3}\right) \cdot \mathbf{E}}{\|\mathbf{E}\|} \frac{\mathbf{N}_{1}}{\left\|\mathbf{N}_{1}\right\|^{2}}-\frac{\left(\mathbf{x}_{2}-\mathbf{x}_{3}\right) \cdot \mathbf{E}}{\|\mathbf{E}\|} \frac{\mathbf{N}_{2}}{\left\|\mathbf{N}_{2}\right\|^{2}}$

$\vb{N}_1 = (\vb{x}_1 - \vb{x}_3)\times(\vb{x}_1 - \vb{x}_4)$ $\vb{N}_2 = (\vb{x}_2 - \vb{x}_4)\times(\vb{x}_2-\vb{x}_3)$ (Cross prod tells the dir. (for normal should be normalized) )

${E} = \vb{x}_4 - \vb{x}_3$

Force:

$\mathbf{f}_{\mathbf{i}}=k \frac{\|\mathbf{E}\|^{2}}{\left\|\mathbf{N}_{1}\right\|+\left\|\mathbf{N}_{2}\right\|} \sin \left(\frac{\pi-\theta}{2}\right) \mathbf{u}_{i}$ $\|\vb{N}_i\|$ ))
$\mathbf{f}_{\mathrm{i}}=k \frac{\|\mathbf{E}\|^{2}}{\left\|\mathbf{N}_{1}\right\|+\left\|\mathbf{N}_{2}\right\|}\left(\sin \left(\frac{\pi-\theta}{2}\right)-\sin \left(\frac{\pi-\theta_{0}}{2}\right)\right) \mathbf{u}_{i}$ $\theta_0$ $\theta_0 = 180^{\circ}$ )

Explicit Derivative is difficult to complete / No energy

A Quadratic Bending Model

2 assumptions: planar case (OK for cloth); little stretching (mainly bending)

Energy function^T $\vb{I}$ is 3x3 identity)

\begin{matrix} E (x) = \frac{1}{2} [\begin{array}{llll} x_{0} & x_{1} & x_{2} & x_{3} \end{array}] Q [\begin{array}{l} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \end{array}]; Q = \frac{3}{A_{0} + A_{1}} q q^{T} \in R^{12 \times 12}; q = [\begin{matrix} (\cot θ_{1} + \cot θ_{3}) I \\ (\cot θ_{0} + \cot θ_{2}) I \\ (- \cot θ_{0} - \cot θ_{1}) I \\ (- \cot θ_{2} - \cot θ_{3}) I \end{matrix}] \in R^{12 \times 3} \\ \Rightarrow E (x) = \frac{3 ∥ q^{T} x ∥^{2}}{2 (A_{0} + A_{1})}; and E (x) = 0 when triangles are flat (co-planar) \end{matrix}

$E(\vb{x}) = 0$

$\vb{Q}$ is a constant matrix => The function is quadratic

Pros & Cons

Easy to implement (even in implicit method)
$\begin{matrix} f (x) = - \nabla E (x) = - Q [\begin{matrix} x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \end{matrix}] H (x) = \frac{\partial^{2} E (x)}{\partial x^{2}} = Q \end{matrix}$
No longer valid if stretches much; Not suitable if the rest config is not planar -> (projective dynamics model / cubic shell model)

The Locking Issue

In the mass-spring / other bending models, assuming cloth planar deformation and cloth bending deformation are independent

In zero bending case: LHS - OK; RHS - For stiff spring NO (Locking Issue) => short of DoFs.

For a manifold mesh, Euler formula: #edges = 3#vertices - 3 - #boundary_edges

If edges are all hard constraints: DoFs = 3 + #boundary_edges

=> no perfect solution

Shape Matching

Lecture 6 Constraint Approaches: PBD / PD / …

Strain Limiting and Position Based Dynamics

The Stiffness Issue

Real-world fabric resist strongly to stretching

But in simulation, only increasing the stiffness can cause: (More expensive / time-costing)

Explicit integrators - Unstable
-> Smaller timesteps and more computational time
Linear systems invloved in implicit integrators will be ill-conditioned
-> More iter and computational time

A Single Spring

If a spring is infinitely stiff -> Length = const

$\phi (\vb{x}) = \|\vb{x}_i - \vb{x}_j\| -L = 0$ $L$ )

$\R^6$ $\vb{x}$ $\delta \boldsymbol{\Omega}$ ) with a shortest path: want the projection

\begin{matrix} {x_{i}^{new}, x_{j}^{new}} = argmin \frac{1}{2} {m_{i} {‖ x_{i}^{new} - x_{i} ‖}^{2} + m_{j} {‖ x_{j}^{new} - x_{j} ‖}^{2}} \\ ⇓ such that ϕ (x) = 0 \\ {\begin{aligned} x^{new} & \leftarrow Projection (x) \\ x_{i}^{new} & \leftarrow x_{i} - \frac{m_{j}}{m_{i} + m_{j}} (‖ x_{i} - x_{j} ‖ - L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖} \\ x_{j}^{new} & \leftarrow x_{j} + \frac{m_{i}}{m_{i} + m_{j}} (‖ x_{i} - x_{j} ‖ - L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖} \end{aligned} \\ ϕ (x^{new}) = ‖ x_{i}^{new} - x_{j}^{new} ‖ - L = ‖ x_{i} - x_{j} - x_{i} + x_{j} + L ‖ - L = 0 \end{matrix}

The 2 new points’ substraction equals the original length => satisfying the constraint

$m_i = m_j$ $m_i = \infty$ for stationary nodes (can be just ignored)

Multiple Springs

Gauss-Seidel Approach

Approach every spring sequentially in a certain order (needs a lot of iter to converge)

Cannot ensure the satisafction of every constraint. More iter give much closer results (to the constraint)
More similar to the stochastic gradient descent (in ML, with spec order)
The order mattrers => cause bias and affect convergence

Jocabi Approach

Projects all edges simultaneously (good for parallization) and then linearly blend the results

sum up all and update together with a weighted average.

Lower convergence rate
More iter give better results

Position Based Dynamics (PBD)

Based on the proj func, similar to shape matching

$\vb{x}$ $\vb{v}$ update as simple particle system (similar to Shape matching in rigid body collision)
PBD: Add constraints
- $\vb{x}$ : use Gauss-Seidel / Jocobi /… to find the projection
- $\vb{v}$ $\vb{x}$ to find

Properties:

The stiffness behavior is subject to non-physical factors
- Iteration num (More iter, slower convergence, less stiffer)
- Mesh resolution (Fewer vertices, faster convergence, stiffer)
The velocity update following projection is important to dynamic effects
Not only use in springs, but triangles, volumes, …

Pros & Cons:

(Usually in 3D softwares)

Pros:
- Parallelable on GPUs (PhysX from NVIDIA)
- Easy to implement
- Fast in low res (Less than 1000 vertices)
- Generic, can handle other coupling and constraints (including fluids)
Cons:
- Not physically correct (no acc sol, stiffness related to meshes account and iter num)
- Low performance in high res
  - Hierarchical approaches (From low to high res. But causes oscillation and other issues)
  - Acceleration approaches (Chebyshev, …)

Strain Limiting

Some normal updates + strain limiting (corrections)

$L$ $\sigma$ - stretching ratio as a limit)

\begin{matrix} x^{new} \leftarrow Projection (x) \\ σ \leftarrow \frac{1}{L} ∥ x_{i} - x_{j} ∥ \\ σ \leftarrow min (max (σ, σ^{min}), σ^{max}) \in [σ^{min}, σ^{max}] \\ {\begin{aligned} x_{i}^{new} & \leftarrow x_{i} - \frac{m_{j}}{m_{i} + m_{j}} (‖ x_{i} - x_{j} ‖ - σ_{0} L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖} \\ x_{j}^{new} & \leftarrow x_{j} - \frac{m_{j}}{m_{i} + m_{j}} (‖ x_{i} - x_{j} ‖ - σ_{0} L) \frac{x_{i} - x_{j}}{‖ x_{i} - x_{j} ‖} \end{aligned} \\ PBD: σ \equiv 1; No limits: σ^{min}, σ^{max} \leftarrow \infty \end{matrix}

Can be used to simulate some cloth whose stiffness increases when being stretched heavily.

The constraints can be used for reducing oscillation for FEM / …

Triangle Area Limit

$s = \sqrt{ \min(\max(A,A^{\min}),A^{\max})/A}$ (Mass center no change)

{x_{i}^{new}, x_{i}^{new}, x_{k}^{new}} = argmin \frac{1}{2} {m_{i} {‖ x_{i}^{new} - x_{i} ‖}^{2} + m_{j} {‖ x_{j}^{new} - x_{j} ‖}^{2} + m_{j} {‖ x_{k}^{new} - x_{k} ‖}^{2}}

Properties:

Avoiding instability and artifacts due to large deformation
Useful for nonlinear effects (Biphasic way)
Help address the locking issue

Projective Dynamics

Main Idea

$\vb{x}$ directly as in PBD, projective dynamics uses projectionquadratic energy $\vb{x}_{e,i}^{\text{new}}$ $\vb{x}_{e,j}^{\text{new}}$ )

\begin{matrix} E (x) = \sum_{e = {i, j}} \frac{k}{2} {‖ (x_{i} - x_{j}) - (x_{e, i}^{new} - x_{e, j}^{new}) ‖}^{2} = \sum_{e = {i, j}} \frac{k}{2} {‖ ∥ x_{i} - x_{j} ∥ \frac{x_{i} - x_{j}}{∥ x_{i} - x_{j} ∥} - L_{e} \frac{x_{i} - x_{j}}{∥ x_{i} - x_{j} ∥} ‖}^{2} = \sum_{e = {i, j}} \frac{k}{2} {(∥ x_{i} - x_{j} ∥ - L_{e})}^{2} \\ {x_{e, i}^{new} - x_{e, j}^{new}} = {Projection}_{e} (x_{i}, x_{j}) for every edge e \end{matrix}

forces $(\vb{x}_{e,i}^{\text{new}} - \vb{x}_{e,j}^{\text{new}}) = \text{const}$ ) (identical to spring force)

f_{i} = - \nabla_{i} E (x) = - k \sum_{e : i \in e} (x_{i} - x_{j}) - (x_{e, i}^{new} - x_{e, j}^{new}) = - k \sum_{e : i \in e} (∥ x_{i} - x_{j} ∥ - L_{e}) \frac{x_{i} - x_{j}}{∥ x_{i} - x_{j} ∥}

=> althought the results are the same, the Hessians are different

\begin{aligned} E (x) = & \frac{k}{2} {‖ (x_{0} - x_{1}) - (x_{0, 0}^{new} - x_{0, 1}^{new}) ‖}^{2} + \frac{k}{2} {‖ (x_{0} - x_{2}) - (x_{1, 0}^{new} - x_{1, 2}^{new}) ‖}^{2} + \frac{k}{2} {‖ (x_{0} - x_{3}) - (x_{2, 0}^{new} - x_{2, 3}^{new}) ‖}^{2} + \\ \frac{k}{2} {‖ (x_{1} - x_{2}) - (x_{3, 1}^{new} - x_{3, 2}^{new}) ‖}^{2} + \frac{k}{2} {‖ (x_{2} - x_{3}) - (x_{4, 2}^{new} - x_{4, 3}^{new}) ‖}^{2} \end{aligned}

f_{0} = [k (x_{0, 0}^{new} - x_{0, 1}^{new}) - k (x_{0} - x_{1})] + [k (x_{1, 0}^{new} - x_{1, 2}^{new}) - k (x_{0} - x_{2})] + [k (x_{2, 0}^{new} - x_{2, 3}^{new}) - k (x_{0} - x_{3})]

\begin{matrix} H = [\begin{matrix} 3 k I & - k I & - k I & - k I \\ - k I & 2 k I & - k I \\ - k I & - k I & 3 k I & - k I \\ - k I & - k I & 2 k I \end{matrix}] \end{matrix}

The Hessian is the 2nd order derivative of energies to the vertices. The coefficients (of the diagonal) can be regarded as the numbers of the springs the vertices connected (3, 2, 3, 3 for vertex 0-3)

$\vb{A = }\frac{1}{\Delta t^2} \vb{M+H}$

=> Applying implicit integration

Simulation by PD

$\vb{A}$ once (usually the most expensive part is the factorization and here the costs can reduce)

Principles

$\vb{H}$ as Hessian

\underset{Hessian}{\underset{⏟}{(\frac{1}{Δ t^{2}} M + H)}} Δ x = \underset{Negative gradient}{\underset{⏟}{- \frac{1}{Δ t^{2}} M (x^{(k)} - x^{[0]} - Δ t v^{[0]}) + f (x^{(k)})}}

$\vb{H}$ usually not calculated exactly (The performance dep. on how well the Hessian is approx., identity matrix can even be used (but convergence takes longer))

Truncating the Hessian to make it p.d. (Lec. 5)
A (magical) diagonal approx (Lab 2)
A const matrix by PD (Lec. 6)

Pros and Cons

Pros

Have theoretical solution with physical meaning
Fast on CPUs with a direct solver (one factorization)
Fast convergence in the first few iter

Cons

Slow on GPUs
Slow convergence over time (fails to consider Hessian caused by projection) (suffering from high stiffness)
Cannot handle constraint changes easily (contacts / remeshing due to fracture / …)

Constrained Dynamics

This method is suitable for strong constraints / contacts / articulated rigid bodies (ragdoll animation) …

A critical problem: what if constraints/forces are very stiff (or even infinitely stiff)
example: multiple rigid bodies: human body can be regarded as several rigid bodies with very strong constraints (joints)

$\phi_e(\vb{x}) = \|\vb{x}_{ei}-\vb{x}_{ej}\| - L_e$

The energy can be rewrited:

\begin{matrix} E (x) = \sum_{e} \frac{1}{2} k {\underset{ϕ_{e}}{\underset{⏟}{(∥ x_{e i} - x_{e j} - L_{e} ∥)}}}^{2} = \frac{1}{2} Φ^{T} (x) C^{- 1} Φ (x) \\ f (x) = - \nabla E = - {(\frac{\partial E}{\partial Φ} \frac{\partial Φ}{\partial x})}^{T} = - J^{T} C^{- 1} Φ = J^{T} λ \end{matrix}

$N$ $E$ $\vb{C}$ $\vb{J}$ $\vb{\lambda}$ - Lagrangian multipliers (Dual variables))

\begin{matrix} Φ (x) \in R^{E}, [\begin{matrix} ϕ_{0} \\ ϕ_{1} \\ ⋮ \\ ϕ_{E} \end{matrix}]; C = [\begin{matrix} \frac{1}{k} \\ \frac{1}{k} \\ ⋱ \end{matrix}] \in R^{E \times E}; J = \frac{\partial Φ}{\partial x} \in R^{E \times 3 N}; λ = - C^{- 1} Φ \in R^{E} \end{matrix}

$\vb{x}$ $\vb{v = \dot {x}}$ $\lambda$ )

\begin{matrix} {\begin{aligned} M v^{new} - Δ t J^{T} λ^{new} = M v \\ C λ^{new} = - Φ^{new} \approx - Φ - J (x^{new} - x) \approx - Φ - Δ t J v^{new} \end{aligned} \\ \Rightarrow [\begin{matrix} M & - Δ t J^{T} \\ Δ t J & C \end{matrix}] [\begin{matrix} v^{new} \\ λ^{new} \end{matrix}] = [\begin{matrix} M v \\ - Φ \end{matrix}] \end{matrix}

Method 1: Solve by a direct solver together (primal-dual fashion): The matrix should be p.d. => need to make p.d.
$\vb{\lambda}^{\text{new}}$ first (not easy to reduce / not always sparse)
$\begin{matrix} (Δ t^{2} {J M}^{- 1} J^{T} + C) λ^{new} = - Φ - Δ t v \\ v^{new} \leftarrow v + - Δ t M^{- 1} J^{T} λ^{new} \end{matrix}$
$\lambda$ first -> get the implicit method formula

$\vb{C}\rightarrow 0$ $\vb{C}$ makes the system better)

Lecture 7 Linear Finite Element Method (FEM)

Linear Finite Element Method (FEM)

The Linear FEM Assumption

$\vb{X}$ $\vb{x = FX + c}$ (“Linear” for this linear formula)

$\vb{F}$ to convert it from ref to deformed state:

\begin{matrix} x_{b a} = x_{b} - x_{a} = F X_{b} + c - F X_{a} - c = F X_{b a} \end{matrix}

For example for edge 10 & 20: (indep on translation)

F X_{10} = x_{10}; F X_{20} = x_{20} \Rightarrow F = [\begin{matrix} x_{10} & x_{20} \end{matrix}] {[\begin{matrix} X_{10} & X_{20} \end{matrix}]}^{- 1}

$\vb{F}$ is related to deformation but contains rotations (want to also reduce rotation)

Green Strain

SVD $\vb{F=UDV}^{\mathrm{T}}$ $\vb{V}^{\mathrm{T}}$ $\vb{D}$ are relevant to deformation.

$\vb{U}$ $\vb{G}$ Green Strain $\vb{F}^{\mathrm{T}}\vb{F} = \vb{VD}\cancel{\vb{U}^{\mathrm{T}}\vb{U}}\vb{DV}^{\mathrm{T}}$ )

\begin{matrix} G = \frac{1}{2} (F^{T} F - I) = \frac{1}{2} (V D^{2} V^{T} - I) = [\begin{matrix} ε_{u u} & ε_{u v} \\ ε_{u v} & ε_{v v} \end{matrix}] \end{matrix}

$\vb{G=0}$ $\|\vb{G}\|$ increases
$\varepsilon_{uu},\varepsilon_{uv}, \varepsilon_{vv}$
$\vb{G}$ $\vb{R}$ $\vb{RF}$ $\vb{G}$ is the same)

Strain Energy Density Function

$W(\vb{G})$ :

Total energy:

E = \int W (G) d A = A^{ref} W (ε_{u u}, ε_{u v}, ε_{v v})

Saint Venant-Kirchhoff Model (StVK) $\lambda$ $\mu$ are Lamé parameters)

W (ε_{u u}, ε_{v v}, ε_{u v}) = \frac{1}{2} {(ε_{u u} + ε_{v v})}^{2} + μ (ε_{u u}^{2} + ε_{v v}^{2} + 2 ε_{u v}^{2})

\begin{matrix} \frac{\partial W}{\partial G} = [\begin{array}{cc} \frac{\partial W}{\partial ε_{u u}} & \frac{1}{2} \frac{\partial W}{\partial ε_{u v}} \\ \frac{1}{2} \frac{\partial W}{\partial ε_{u v}} & \frac{\partial W}{\partial ε_{v v}} \end{array}] = [\begin{array}{cc} 2 μ ε_{u u} + λ ε_{u u} & 2 μ ε_{u v} \\ 2 μ ε_{u v} & 2 μ ε_{v v} + λ ε_{v v} \end{array}] = 2 μ G + λ trace (G) I = S \end{matrix}

$\vb{S}$ - Second Piola-Kirchhoff stress tensor (some form of force density))

Force:

f_{i} = - {(\frac{\partial E}{\partial x_{i}})}^{T} = - A^{ref} {(\frac{\partial W}{\partial x_{i}})}^{T} = - A^{ref} {(\frac{\partial W}{\partial ε_{u u}} \frac{\partial ε_{u u}}{\partial x_{i}} + \frac{\partial W}{\partial ε_{v v}} \frac{\partial ε_{v v}}{\partial x_{i}} + \frac{\partial W}{\partial ε_{u v}} \frac{\partial ε_{u v}}{\partial x_{i}})}^{T}

$\frac{\partial W}{\partial \varepsilon_{ij}}$ are solved before, we concentrate more on the next part

$\mathbf{F}=\left[\begin{array}{ll} \mathbf{x}_{10} & \mathbf{x}_{20} \end{array}\right]\left[\begin{array}{ll} \mathbf{r}_{10} & \mathbf{r}_{20} \end{array}\right]^{-1}=\left[\begin{array}{ll} \mathbf{x}_{10} & \mathbf{x}_{20} \end{array}\right]\left[\begin{array}{ll} a & b \\ c & d \end{array}\right]=\left[\begin{array}{ll} a \mathbf{x}_{10}+c \mathbf{x}_{20} & b \mathbf{x}_{10}+d \mathbf{x}_{20} \end{array}\right]$

By def

\begin{matrix} G = \frac{1}{2} (F^{T} F - I) = [\begin{array}{cc} \frac{1}{2} {(a x_{10} + c x_{20})}^{T} (a x_{10} + c x_{20}) - \frac{1}{2} & \frac{1}{2} {(a x_{10} + c x_{20})}^{T} (b x_{10} + d x_{20}) \\ \frac{1}{2} {(a x_{10} + c x_{20})}^{T} (b x_{10} + d x_{20}) & \frac{1}{2} {(b x_{10} + d x_{20})}^{T} (b x_{10} + d x_{20}) - \frac{1}{2} \end{array}] \end{matrix}

\begin{matrix} \Rightarrow {\begin{aligned} \frac{\partial ε_{u u}}{\partial x_{1}} = a {(a x_{10} + c x_{20})}^{T} \frac{\partial ε_{v v}}{\partial x_{1}} = b {(b x_{10} + d x_{20})}^{T} \frac{\partial ε_{u v}}{\partial x_{1}} = \frac{1}{2} a {(b x_{10} + d x_{20})}^{T} + \frac{1}{2} b {(a x_{10} + c x_{20})}^{T} \\ \frac{\partial ε_{u u}}{\partial x_{2}} = c {(a x_{10} + c x_{20})}^{T} \frac{\partial ε_{v v}}{\partial x_{2}} = d {(b x_{10} + d x_{20})}^{T} \frac{\partial ε_{u v}}{\partial x_{2}} = \frac{1}{2} c {(b x_{10} + d x_{20})}^{T} + \frac{1}{2} d {(a x_{10} + c x_{20})}^{T} \end{aligned} \end{matrix}

(the method is too complex => simplification)

\begin{aligned} f_{1} & = - A^{ref} {(\frac{\partial W}{\partial ε_{u u}} \frac{\partial ε_{u u}}{\partial x_{1}} + \frac{\partial W}{\partial ε_{v v}} \frac{\partial ε_{v v}}{\partial x_{1}} + \frac{\partial W}{\partial ε_{u v}} \frac{\partial ε_{u v}}{\partial x_{1}})}^{T} \\ = - A^{ref} (\frac{\partial W}{\partial ε_{u u}} a {(a x_{10} + c x_{20})}^{T} + \frac{\partial W}{\partial ε_{v v}} b {(b x_{10} + d x_{20})}^{T} + \frac{\partial W}{\partial ε_{u v}} \frac{1}{2} a {(b x_{10} + d x_{20})}^{T} + \frac{\partial W}{\partial ε_{u v}} \frac{1}{2} b {(a x_{10} + c x_{20})}^{T}) \\ = - A^{ref} ([\begin{array}{cc} a x_{10} + c x_{20} & b x_{10} + d x_{20} \end{array}] [\begin{array}{cc} \frac{\partial W}{\partial ε_{u u}} a + \frac{\partial W}{\partial ε_{u v}} \frac{1}{2} b \\ \frac{\partial W}{\partial ε_{u v}} \frac{1}{2} a + \frac{\partial W}{\partial ε_{v v}} b \end{array}]) \\ = - A^{ref} \underset{deformation gradient}{\underset{⏟}{[a x_{10} + c x_{20} b x_{10} + d x_{20}]}} \underset{Second Piola-Kirchhoff Stress}{\underset{⏟}{[\begin{array}{cc} \frac{\partial W}{\partial ε_{u u}} & \frac{1}{2} \frac{\partial W}{\partial ε_{u v}} \\ \frac{1}{2} \frac{\partial W}{\partial ε_{u v}} & \frac{\partial W}{\partial ε_{v v}} \end{array}]}} [\begin{matrix} a \\ b \end{matrix}] = - A^{ref} F S [\begin{array}{c} a \\ b \end{array}] \end{aligned}

\Rightarrow [f_{1} f_{2}] = - A^{ref} F S [r_{10} r_{20}]^{- 1}; f_{0} = - f_{1} - f_{2}

For tetrahedron (3D ref -> 3D deformation)

Finite Volume Method (FVM)

For simple triangle / tetrahedra, FVM is equivalent to linear FEM

FVM considers force calculation in an integration perspective other than a differentiation perspective

Traction and Stress in FVM

Basic Ideas in FVM

General Case:

Traction $\vb{t}$ total force at the interface $\vb{f } = \oint_L \vb{t}\ \mathrm{d}l$
$\vb{\sigma}$ $\vb{t = \sigma n}$ (projection: interface normal -> traction)

The Triangle Case $\vb{x}_0$ )

Force contributed by an element $\vb{f}_0 = \oint _L \vb{\sigma n }\ \mathrm{d}l$ (the integral of traction)
Divergence theorem $\oint_L \vb{\sigma n}\ \mathrm{d}l + \oint_{L_{20}} \vb{\sigma n}\ \mathrm{d}l + \oint_{L_{10}}\vb{\sigma n}\ \mathrm{d}l = 0$ $\vb{\sigma}$ is constant within the elem)
Force $\vb{f}_0 = -\oint_{L_{20}}\vb{\sigma n}_{20}\ \mathrm{d}l -\oint_{L_{10}}\vb{\sigma n}_{10}\ \mathrm{d}l = -\vb{\sigma} (\frac{\|\vb{x}_{20}\|}{2}\vb{n}_{20} + \frac{\|\vb{x}_{10}\|}{2}\vb{n}_{10})$ $\sigma$ is unknown)

In 3D (tetrahedra): FVM does the surfacial integral

Force:
$\begin{aligned} f_{0} & = - \oint_{Ω}^{σ} n d A = - σ (\frac{A_{012}}{3} n_{012} + \frac{A_{023}}{3} n_{023} + \frac{A_{031}}{3} n_{031}) \\ = - \frac{σ}{3} (\frac{‖ x_{10} \times x_{20} ‖}{2} \frac{x_{10} \times x_{20}}{‖ x_{10} \times x_{20} ‖} + \frac{‖ x_{20} \times x_{30} ‖}{2} \frac{x_{20} \times x_{30}}{‖ x_{20} \times x_{30} ‖} + \frac{‖ x_{30} \times x_{10} ‖}{2} \frac{x_{30} \times x_{10}}{‖ x_{30} \times x_{10} ‖}) \\ = - \frac{σ}{6} (x_{10} \times x_{20} + x_{20} \times x_{30} + x_{30} \times x_{10}) \end{aligned}$

Different Stresses

The stress tensor here: mapping from the interface normal to the traction (but kind of different from FEM)

$\vb{S}$ $\vb{N}$ $\vb{T}$ are in the reference state
$\vb{n}$ $\vb{t}$ are in the deformed state

Output \ Input	$\vb{N}$ in the ref state (Undeformed)	$\vb{n}$ in the current state (Deformed)
*Traction in the ref* state (unformed)**	$\vb{S}$ )	-
*Traction in the cur* state (formed)**	$\vb{P}\vb{= FS}$ )	$\vb{\sigma} = \det ^{-1} \vb{(F)PF}^{\mathrm{T}} = \det^{-1}\vb{(F)FSF^{\mathrm{T}}}$ )

$\vb{F}$ - Deformation gradient)

The Finite Volume Method

\begin{aligned} f_{1} & = - \frac{σ}{6} (x_{01} \times x_{21} + x_{21} \times x_{31} + x_{31} \times x_{01}) (Replace σ with P (=> normal)) \\ = - \frac{P}{6} (X_{01} \times X_{21} + X_{21} \times X_{31} + X_{31} \times X_{01}) (The term including X is constant => precomputed) \\ = - \frac{F S}{6} b_{1} (b_{1} is the constant) \end{aligned}

$\vb{S} = \frac{\partial W}{\partial \vb{G}}$ as in previous FEM formulation

$\vb{b}_1$ $\vb{X}_1$ $\vb{b}_1$ )

$\vb{X}_1$

\begin{aligned} X_{10}^{T} b_{1} & = X_{10}^{T} (X_{01} \times X_{21} + X_{21} \times X_{31} + X_{31} \times X_{01}) = X_{10}^{T} (X_{21} \times X_{31}) \\ = X_{01}^{T} (X_{31} \times X_{21}) = 6 Vol \end{aligned}

$\vb{X}_{20}^{\mathrm{T}}$ $\vb{X}_{30}^{\mathrm{T}}$ equal to 0

\begin{matrix} \Rightarrow {[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{T} b_{1} = [\begin{matrix} 6 Vol \\ 0 \\ 0 \end{matrix}]; {[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{T} b_{2} = [\begin{matrix} 0 \\ 6 Vol \\ 0 \end{matrix}]; {[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{T} b_{3} = [\begin{matrix} 0 \\ 0 \\ 6 Vol \end{matrix}] \\ \Rightarrow [\begin{matrix} b_{1} & b_{2} & b_{3} \end{matrix}] = 6 Vol {[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{- T} = \frac{1}{\det ({[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{- 1})} {[\begin{matrix} X_{10} & X_{20} & X_{30} \end{matrix}]}^{- T} \end{matrix}

Algorithm (Explicit)

-> Various material (Heterogeneous materials) / various elastic / …

Hyperelastic Models

StVK: cannot handle flips / … (only widely used in graphics)
Hyperelastic: More used in engineering (material / mechanics /… ) => more general

Isotropic Materials

Isotropic

Force on any direction affects the same (with same behaviors) (-> ansotropic)

$\vb{F}$ $\vb{D} = \lambda_0, \lambda_ 1, \lambda_2$ Principal stretches $\vb{F}$

P (F) = P (U D V^{T}) = U P (λ_{0}, λ_{1}, λ_{2}) V^{T}

$\vb{P}(I_{\vb{C}},II_{\vb{C}}, III_{\vb{C}})$ $\vb{C} = \vb{F}^{\mathrm{F}}\vb{U}$ is the right Cauchy-Green deformation tensor)

\begin{matrix} I_{C} = tr (C) = λ_{0}^{2} + λ_{1}^{2} + λ_{2}^{2}; I I I_{C} = tr (C^{2}) = λ_{0}^{4} + λ_{1}^{4} + λ_{2}^{4} \\ I I_{C} = \frac{1}{2} ({tr}^{2} (C) - tr (C^{2})) = λ_{0}^{2} λ_{1}^{2} + λ_{0}^{2} λ_{2}^{2} + λ_{1}^{2} λ_{2}^{2} \end{matrix}

Models

The Saint Venant-Kirchhoff model (StVK):
$W = \frac{s_{0}}{2} {(I_{C} - 3)}^{2} + \frac{s_{1}}{4} (I I_{C} - 2 I_{C} + 3)$
The neo-Hookean model: (More standard model in material mechanics)
$W = s_{0} (I I I_{C}^{- 1 / 3} I_{C} - 3) + s_{1} (I I I_{C}^{- 1 / 2} - 1)$
$W$ : left - Against shearing; right - Against bulky change (volume change) (same as StVK)
The Mooney Rivlin model:
$W = s_{0} (I I I_{C}^{- 1 / 3} I_{C} - 3) + s_{1} (I I I_{C}^{- 1 / 2} - 1) + s_{2} (\frac{1}{2} I I I_{C}^{- 2 / 3} (I_{C}^{2} - I I_{C}) - 3)$
The Fung model:
$W = s_{0} (I I I_{C}^{- 1 / 3} I_{C} - 3) + s_{1} (I I I_{C}^{- 1 / 2} - 1) + s_{2} (e^{s_{3} (I I I_{C}^{- 1 / 3} I_{C} - 3)} - 1)$

$\vb{P}$ )

Use the principal stretch for computation => first P-K stress

\begin{matrix} P (λ_{0}, λ_{1}, λ_{2}) = [\begin{matrix} \frac{\partial W}{\partial λ_{0}} \\ \frac{\partial W}{\partial λ_{1}} \\ \frac{\partial W}{\partial λ_{2}} \end{matrix}] \Rightarrow P = U P (λ_{0}, λ_{1}, λ_{2}) V^{T} \end{matrix}

Algorithm

$\vb{P} = \vb{F}\frac{\partial W}{\partial \vb{G}}$

The Limitation of StVK

The compression (with high stress) and flip cannot be handled with StVK model (even get stable in the flipped side) => UNSTABLE

=> Neo Hookean solve the problem

Lecture 8 Linear Finite Element Method II

Nonlinear Optimization

Gradient Descent

$\vb{x}^* = \operatorname{argmin} F(\vb{x})$ => The negative gradient dir is the fastest descent dir

Algorithm

Step Size Adjustment

$\beta$ is a smaller num (~ 0.3-0.4) -> Wolfe Condition)

Descent Directions

$\vb{d(x)}$ $\alpha$ $F(\vb{x}) > F(\vb{x+ \alpha d(x)} )$ (at least descending)

$-\grad F(\vb{x})\cdot \vb{d(x)}>0$ (the same side with the negative dir)

$F(\vb{x}^{(0)}) > F(\vb{x}^{(1)}) > F(\vb{x}^{(2)}) > F(\vb{x}^{(3)}) > \dots$

Descent Methods

Gradient descent $\vb{d(x)} = -\grad F(\vb{x}) \Rightarrow -\grad F(\vb{x}) \cdot (-\grad F(\vb{x})) >0$
Newton’s method (if the Hessian is always positive definite):
$d (x) = - {(\frac{\partial^{2} F (x)}{\partial x^{2}})}^{- 1} \nabla F (x) \Rightarrow - \nabla F (x) \cdot (- {(\frac{\partial^{2} F (x)}{\partial x^{2}})}^{- 1} \nabla F (x)) > 0$
positive definite $\vb{P}$ to modify the gradient yields a descent method:
$d (x) = - P^{- 1} \nabla F (x) \Rightarrow - \nabla F (x) \cdot (- P^{- 1} \nabla F (x)) > 0$

A unified descent framework:

Differences between these methods:

Total cost = Per-iter cost * Num of iter

Lecture 9 Collision Handling

Collision Detection

Collision Detection Pipeline

Broad-Phase: Remove the pairs which are impossible to collide
- SH: Easy to implement; GPU friendly; needs to recompute after updating obj
- BVH: More involved; not GPU friendly; update BV to update BVH
Narrow-Phase: Which pairs really collide

Board-Phase Collision Culling

Spatial Hashing (Partitioning)

Static space division

Spatial partitioning divides the space by a grid and stores objects into grid cells

$t_4$	$t_4$	$t_5$	$t_5$
$t_2$	$t_3$	$t_3$ $t_5$	Cell 11: N/A
$t_0$ $t_2$	$t_0$ $t_3$	$t_3$	$t_1$
$t_0$	Cell 1: N/A	$t_1$	$t_1$

$t_3$ $t_0$ $t_5$ $t_0t_3$ $t_3t_5$ and do further detection

For moving objects, just expand the object region (add more)

Optimization for memory costs: Memory wastes in empty cells and too many meshes should be built in 3D. Instead of allocating memories to cell, build an object-cell list and sort them

e.g., (not relavent to space)


xxxxxxxxxx
{0, t_0}; {4, t_0}; {5, t_0}; {2, t_1}; {3, t_1}; {7, t_1}; {4, t_2}; {8, t_2}; {5, t_3}; {6, t_3}; {9, t_3}; {10, t_3}; {12, t_4};  {13, t_4}; {10, t_5}; {14, t_5}; {15, t_5}

Then sort by cell ID => easy to find out what is stored in each cell


xxxxxxxxxx
{0, t_0}; {2, t_1}; {3, t_1}; {4, t_0}; {4, t_2}; {5, t_0}; {5, t_3}; {6, t_3}; {7, t_1}; {8, t_2}; {9, t_3}; {10, t_3}; {10, t_5}; {12, t_4};  {13, t_4}; {14, t_5}; {15, t_5}

Morton Code

One question is how to define the cell ID. Using the grid order is not optimal, since it cannot be easily extended and it is lack of locality (LHS). Morton code uses a Z-pattern instead (RHS).

Bounding Volume Hierarchy (BVH)

Bounding volume hierarchy is built on geometric/topological proximity of objects.

External Object

The parts that are not like to collide in different boxes => test if box intersection occurs

To find elem potentially in collision with an obj, just traverse the tree:

Internal Intersections (Self Collision)

2 procedures in the tree


xxxxxxxxxx
Process_Node(A)
{   // check if the parental node (A) has self-collision
    For every A’s child: B
        Process_Node(B)
    For every A’s children pair <B, C>
        if B and C intersect
            Process_Pair(B, C)  // the next precedure
}


xxxxxxxxxx
Process_Pair(B, C)
{   // check if the children pair intersects
    For every B’s child: B’
    For every C’s child: C’
        if B’ and C’ intersect
            Process_Pair(B’, C’)
}

BVH Models

Axis-aligned bounding box (AABB) / Spheres / Oriented bounding box (OBB)

Two AABBs intersect if and only if they intersect in every axis; shperes: need to compute the distance

Narrow-Phase Collision Test

Discrete Collision Detection (DCD)

Actually not collision detection but intersection detection

discrete time instant $\vb{x}^{[0]}, \vb{x}^{[1]}, \dots$

To a triangle mesh, the basic test is edge-triangle intersection test: easy and robust

$t$ means the interpolation value between the points other than time)

Tunneling

The difference between collision and intersection detection:

Intersection detection only considers if the 2 triangles (/…) have intersection at the two discrete time steps
Collision detection can consider the penentration

Problem $\vb{x}^{[0]}$ $\vb{x}^{[1]}$ doesn’t have intersection with the blue one but actually the position relationship changes => tunneling problem (penetrating through each other)

Not usually occurs for bulky objects but can happen for thin surfaces

Continuous Collision Detection (CCD)

between two states $\vb{x}^{[0]}$ $\vb{x}^{[1]}$

To a triangle mesh, there two basic tests: vertex-triangle and edge-edge tests

Vertex-Triangle:

$\vb{x}_{30} = \vb{x}_3 - \vb{x}_{0}$ which is a function about time)

Need to solve a third order equation (with one variable) => has formulae solution (not recommended, need to solve cubic root) / use binary division (0, 1) and get the first sol

Edge-Edge:

Issues with CCD:

Floating-point errors, especially due to root finding of a cubic equation
- Buffering epsilons, but that causes false positives
- Gaming GPUs often use single floating-point precision
Computational costs: more expensive than DCD
- Some argue that broad-phase collision culling is the bottleneck
Difficulty in implementation

Continuous Collision Response Approaches

The Two Continuous Collision Response Approaches

$\vb{x}^{[1]}$ $\overline{\vb{x}}^{[1]}$ $\vb{x}^{[0]}$ $\overline {\vb{x}}^{[1]}$ is intersection-free

Interior Point Method: Use the rest state to the collided state to estimate the every next state (safer)
- Slow: far from solution; all vertices; cautiously by small step size
- Always succeed
Impact Zone Method: Use the final point to estimate, then do optimization in the “impact zone”
- Fast: Close to solution; only vertices in collision (impact zones); can take large step sizes
- May not succeed

Interior Point Methods (Cont.)

Log-Barrier Interior Point Methods

For simplicity, just consider the Log-barrier repulsion between 2 vertices (usually in graphics do a cutoff at some distance)

\begin{matrix} E (x) = - ρ \log ∥ x_{i j} ∥ \\ f_{i} (x) = - \nabla_{i} E = ρ \frac{x_{i j}}{∥ x_{i j} ∥^{2}}; f_{j} (x) = - \nabla_{j} E = - ρ \frac{x_{i j}}{∥ x_{i j} ∥^{2}} \end{matrix}

Implementation

Formulate the problem:

{\overset{―}{x}}^{[1]} \leftarrow {argmin}_{x} (\frac{1}{2} ∥ x - x^{[1]} ∥^{2} - ρ \sum \log ∥ x_{i j} ∥)

Gradient Descent:

$\alpha$ $\alpha$ , we need collision tests.

Impact Zone Optimization (Cont.)

$\vb{x}^{[1]}$ until it becomes intersection-free. (This potentially suffers from the tunneling issue, uncommon)

{\overset{―}{x}}^{[1]} \leftarrow {argmin}_{x} \frac{1}{2} ∥ \overset{―}{x} - x^{[1]} ∥^{2}

Such that:

{\begin{cases} C (x) = - (x_{3} - b_{0} x_{0} - b_{1} x_{1} - b_{2} x_{1}) \cdot N \leq 0 & For each detected vertex-triangle pair \\ C (x) = - (b_{2} x_{2} + b_{3} x_{3} - b_{0} x_{0} - b_{1} x_{1}) \cdot N \leq 0 & For each detected edge-edge pair \end{cases}

The optimizations are always trying to get close to the solution

Rigid Impact Zones

The rigid impact zone method simply freezes vertices in collision from moving in their pre-collision state. It’s simple and safe, but has noticeable artifacts

Only used when all other methods not work

A Practical System

Untangling Cloth (Discrete)

Intersection Elimination

Consider how to eliminate existing intersections, but without using any collision history

Useful when there are already intersections in simulation, due to:

Past collision handling failures
Intense user interaction

In this case, we don’t require the simulation is to always intersection-free

Eliminating cloth-volume and volume-volume intersections is straightforward: simply pushing vertices/edges in the volume out

Untangling Cloth

Cloth-cloth intersection is complicated (don’t have a clear definition of inside and outside)

Baraff et al. used flood-fill to segment cloth into regions and decided which region is in intersection. (Cannot handle boundary well, not friendly in GPU)
divide into several regions => which is larger => the smaller part moves to the larger side
Volino and Magnenat-Thalmann proposed to untangle cloth by reducing the intersection contour
Their method can handle boundaries, but it doesn’t always work

Lecture 10 Surface Waves

Two Types of Simulation Approaches

Lagrangian Approach: Dynamic particles or mesh (Node movement carries physical quantities (mass, velocity, …))
Eulerian Approach: Static grid or mesh (Grid/Mesh doesn’t move; Stored physical quantities change)

A Height Field Model

Height Field

$h(x)$ $h$ $u(x)$ (negative -> opposite direction)

$h\cdot u$ means the volume flowing through a position at the unit time)

\begin{matrix} {\begin{aligned} \frac{d h (x)}{d t} + \frac{d (h (x) u (x))}{d x} = 0 \\ \frac{d u (x)}{d t} = \underset{advection}{\underset{⏟}{- u (x) \frac{d u (x)}{d x}}} - \frac{1}{ρ} \frac{d P (x)}{d x} + \underset{ext.}{\underset{⏟}{a (x)}} \end{aligned} \end{matrix}

Ignoring advectionexternal acceleration $\rho$ $P(x)$ - pressure)

\frac{d u (x)}{d t} = - \frac{1}{ρ} \frac{d P (x)}{d x}

$x$ $x+\mathrm{d}x$ $\mathrm{d}P(x) = P(x+\mathrm{d}x) - P(x)$ => the velocity changes due to the change of pressure (also dep on density)

$u(x)$ $\mathrm{d}P(x)<0$ $\mathrm{d}P(x)>0$

Shallow Wave Equation

Simplification of the two equations by assuming shallow wave (the wave is small and the height change can be neglected)

\begin{matrix} {\begin{aligned} \frac{d h}{d t} + \frac{d (h u)}{d x} = 0 ⟶ \frac{d h}{d t} + u \frac{d h}{d x} + h \frac{d u}{d x} = 0 \\ \frac{d u}{d t} = - \frac{1}{ρ} \frac{d P}{d x} \end{aligned} \begin{aligned} \overset{d t}{⟶} \\ \overset{d x}{⟶} \end{aligned} {\begin{aligned} \frac{d^{2} h}{d t^{2}} + h \frac{d^{2}}{d x d t} = 0 \\ \frac{d^{2} u}{d x d t} = - \frac{1}{ρ} \frac{d^{2} P}{d x^{2}} \end{aligned} \end{matrix}

$u$ and formulate the shallow wave equation (indep on the velocity field)

\Rightarrow \frac{d^{2} h}{d t^{2}} = \frac{h}{ρ} \frac{d^{2} P}{d x^{2}}

Finite Differencing

Discretization

Finite Difference

Use the difference to approximate the derivative

Forward differencing (first order, 2nd order error)
$f (t_{0} + Δ t) = f (t_{0}) + Δ t \frac{d f (t_{0})}{d t} + \frac{Δ t^{2} d^{2} f (t_{0})}{2} + \dots \Rightarrow \frac{d f (t_{0})}{d t} \approx \frac{f (t_{0} + Δ t) - f (t_{0})}{Δ t}$
Backward differencing (first order)
$f (t_{0} - Δ t) = f (t_{0}) - Δ t \frac{d f (t_{0})}{d t} + \frac{Δ t^{2} d^{2} f (t_{0})}{2} + \dots \Rightarrow \frac{d f (t_{0})}{d t} \approx \frac{f (t_{0}) - f (t_{0} - Δ t)}{Δ t}$
Central differencing (second order)
$\Rightarrow \frac{d f (t_{0})}{d t} \approx \frac{f (t_{0} + Δ t) - f (t_{0} - Δ t)}{2 Δ t}$

Second-Order Derivatives (central differencing in the fluid field)

height $\mathrm{d}^2h_i / \mathrm{d}t^2$ (Laplacian)

\frac{d^{2} h_{i} (t_{0})}{d t^{2}} \approx \frac{\frac{d h_{i} (t_{0} + 0.5 Δ t)}{d t} - \frac{d h_{i} (t_{0} - 0.5 Δ t)}{d t}}{Δ t} \approx \frac{h_{i} (t_{0} + Δ t) + h_{i} (t_{0} - Δ t) - 2 h_{i} (t_{0})}{Δ t^{2}}

pressure $\mathrm{d}^2P / \mathrm{d}x^2$

\frac{d^{2} P_{i}}{d x^{2}} \approx \frac{\frac{d P_{i + 0.5}}{d x} - \frac{d P_{i - 0.5}}{d x}}{Δ x} \approx \frac{P_{i + 1} + P_{i - 1} - 2 P_{i}}{Δ x^{2}}

Discretized Shallow Wave Equation

Want to obtain the height field in the next time step

\frac{d^{2} h}{d t^{2}} = \frac{h}{ρ} \frac{d^{2} P}{d x^{2}} \Rightarrow h_{i} (t_{0} + Δ t) = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i})

Problem: Volume loss or generation

Volume Preservation

We want the volume to stay the same

$\sum h_i(t) = \sum h_i(t-\Delta t) = V$ , but by using the formuation in the prev section, we will get

\begin{aligned} \sum h_{i} (t + Δ t) & = 2 {\sum h_{i} (t_{0})}^{V} - {\sum h_{i} (t_{0} - t)}^{V} + \sum \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \\ = V + \sum \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \end{aligned}

$\sum \frac{\Delta t^2 h_i}{\Delta x^2 \rho} (P_{i+1} + P_{i-1} -2P_i)$ is 0

Solution 1

$h_i$ $\frac{h_{i-1}+h_i}{2}$ $\frac{h_{i+1}+h_i}{2}$ instead)

\begin{aligned} h_{i} (t_{0} + Δ t) & = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} h_{i}}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \\ \Rightarrow h_{i} (t_{0} + Δ t) & = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2}}{Δ x^{2} ρ} ((\frac{h_{i - 1} + h_{i}}{2}) (P_{i - 1} - P_{i}) + (\frac{h_{i + 1} + h_{i}}{2}) (P_{i + 1} - P_{i})) \end{aligned}

$h_i$ $h_{i+1}$ . The water flowing from the LHS will flow into the RHS => Volume preserved

Solution 2

$h_i$ in the right term is const

\begin{aligned} \Rightarrow h_{i} (t_{0} + Δ t) & = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \frac{Δ t^{2} H}{Δ x^{2} ρ} (P_{i + 1} + P_{i - 1} - 2 P_{i}) \\ So that \sum h_{i} (t_{0} + Δ t) & = V + \frac{Δ t^{2} H}{Δ x^{2} ρ} \sum \underset{Must be 0}{\underset{⏟}{((P_{i + 1} - P_{i}) + (P_{i - 1} - P_{i}))}} \end{aligned}

Pressure

height $P_i = \rho g h_i$

\Rightarrow h_{i} (t_{0} + Δ t) = 2 h_{i} (t_{0}) - h_{i} (t_{0} - Δ t) + \underset{α}{\underset{⏟}{\frac{Δ t^{2} H g}{Δ x^{2}}}} (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0}))

$\alpha$

Viscosity

$\beta$ : too high: not viscos; too low: not like fluid)

\begin{matrix} h_{i} (t_{0} + Δ t) = h_{i} (t_{0}) + β \underset{Momentum here}{\underset{⏟}{(h_{i} (t_{0}) - h_{i} (t_{0} - Δ t))}} + α (h_{i + 1} (t_{0}) + h_{i - 1} (t_{0}) - 2 h_{i} (t_{0})) \end{matrix}

Boundary Conditions

Dirichlet $H_{i+1}$ is constant (open boundary)
$h_{i + 1} - h_{i} + h_{i - 1} - h_{i} = H_{i + 1} - h_{i} + h_{i - 1} - h_{i}$
Neumann Boundary: specifies the boundary derivativeszero-derivative $h_{i+1}\equiv h_i$ (closed boundary)
$h_{i + 1} - h_{i} + h_{i - 1} - h_{i} = h_{i - 1} - h_{i}$

Algorithm with Neumann Boundary

– to 3D –

Two-way Coupling

The coupling will be two ways: between liquid->solid and solid->liquid (also can be bubbles / soft solids / …)

Key: how to expel water out of the gray cell regions

Virtual Height

Idea $v_i$ $h_i^{\text{real\_new}} = h_i - e_i$

\begin{matrix} {\begin{aligned} Left: & h_{i} - e_{i} = h_{i} + β (h_{i} - h_{i}^{old}) + α (ν_{i + 1} + h_{i + 1} + h_{i - 1} - 2 ν_{i} - 2 h_{i}) = h_{i}^{new} + α (ν_{i + 1} - 2 ν_{i}) \\ Right: & h_{i + 1} - e_{i + 1} = h_{i + 1} + β (h_{i + 1} - h_{i + 1}^{old}) + α (h_{i + 2} + ν_{i} + h_{i} - 2 ν_{i + 1} - 2 h_{i + 1}) = h_{i + 1}^{new} + α (ν_{i} - 2 ν_{i + 1}) \end{aligned} \end{matrix}

Poisson’s Equation

$\nu_i$ $\nu_{i+1}$ being unknowns

\begin{aligned} 2 ν_{i} - ν_{i + 1} & = \frac{1}{α} (h_{i}^{new} - h_{i} + e_{i}) = b_{i} \\ - ν_{i} + 2 ν_{i + 1} & = \frac{1}{α} (h_{i + 1}^{new} - h_{i + 1} + e_{i + 1}) = b_{i + 1} \end{aligned}

$\nu_{i-1}$ $\nu_{i+2}$ are all 0, no virtual height)

[\begin{matrix} 2 & - 1 \\ - 1 & 2 \end{matrix}] [\begin{matrix} ν_{i} \\ ν_{i + 1} \end{matrix}] = [\begin{matrix} b_{i} \\ b_{i + 1} \end{matrix}] \Rightarrow [\begin{array}{cccc} 1 \\ - 1 & 2 & - 1 \\ - 1 & 2 & - 1 \\ 1 & 1 \end{array}] = [\begin{matrix} ν_{i - 1} \\ ν_{i} \\ ν_{i + 1} \\ ν_{i + 2} \end{matrix}] = [\begin{matrix} 0 \\ b_{i} \\ b_{i + 1} \\ 0 \end{matrix}]

Algorithm

Process: Set boundary condition -> MaskPCG_Solve() $\nu$ -> Update the height field (Usually need to multiply a factor, prevent too large waves when dragging the cube <- explicit’s unstability)

Rigid Body Update

Estimate the floating (buoyancy) force by the actual water expelled in every column

\begin{matrix} f_{i} = ρ g Δ x (h_{i} - h_{i}^{new}) & in 2D \\ f_{i, j} = ρ g Δ A (h_{i, j} - h_{i, j}^{new}) & in 3D \end{matrix}

Also need to consider rotation => torque

Lecture 11 Eulerian Fluids

A Grid Representation and Finite Differencing

A Regular Grid Representation

Storing physical quantities over a grid => Fields

Scalars: Density/color, Pressure, Temperature, …
Vectors: Velocities

Finite Differencing on Grid

Grid is good for central differencing

\begin{matrix} \begin{array}{r} \frac{\partial f_{i + 0.5, j}}{\partial x} \approx \frac{f_{i + 1, j} - f_{i, j}}{h} \\ \frac{\partial f_{i - 0.5, j}}{\partial x} \approx \frac{f_{i, j} - f_{i - 1, j}}{h} \end{array}} \Rightarrow \frac{\partial^{2} f_{i, j}}{\partial x^{2}} \approx \frac{\frac{\partial f_{i + 0.5, j}}{\partial x} - \frac{\partial f_{i - 0.5, j}}{\partial x}}{h} \approx \frac{f_{i - 1, j} + f_{i + 1, j} - 2 f_{i, j}}{h^{2}} \end{matrix}

$y$ direction)

Discretized Laplacian

Δ f_{i, j} = \frac{\partial^{2} f_{i, j}}{\partial x^{2}} + \frac{\partial^{2} f_{i, j}}{\partial y^{2}} \approx \frac{f_{i - 1, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1} - 4 f_{i, j}}{h^{2}}

Boundary Conditions

$f_{i-1,j}$ if it’s outside

Dirichlet $f_{i-1,j} = C$
$Δ f_{i, j} \approx \frac{C + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1} - 4 f_{i, j}}{h^{2}}$
Neumann $f_{i-1,j} =f_{i,j}$
$Δ f_{i, j} \approx \frac{f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1} - 3 f_{i, j}}{h^{2}}$

Example: Laplace’s Equation

$\Delta f = 0$ $f$ , every )

f_{i - 1, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1} - 4 f_{i, j} = 0

Note: at least one condition should be Dirichlet condition (avoid inifinite solutions)

Laplacian Smoothing is basically average the self value with neighbors => the distribution will be smoother and smoother

$\alpha$ term)

The process of applying Laplacian smoothing is called diffusion

Staggered Grid

Problem with Central Differencing

Central differencing gives the derivative in the middle

$(i+0.5, j)$
$\frac{\partial f_{i,j}}{\partial x}$ $f_{i-1,j}$ $f_{i+1,j}$ $f_{i,j}$ is unused

Solution: Staggered Grid

We define some physical quantities on faces, specifically velocities

The x-part of the velocity is defined on vertical faces
The y-part of the velocity is defined on horizonal faces

flow speed between two cells $(i,j)$ $u_{i+1,j} + v_{i,j+1} - u_{i,j} - v_{i,j}$

Divergence-Free Condition

No volume change is equal to say the fluid is incompressible (divergence-free velocity field)

u_{i + 1, j} + v_{i, j + 1} - u_{i, j} - v_{i, j} = 0 \Rightarrow \nabla \cdot u_{i, j} = \frac{\partial u_{i, j}}{\partial x} + \frac{\partial v_{i, j}}{\partial y} \approx \frac{u_{i + 1, j} - u_{i, j}}{h} + \frac{v_{i, j + 1} - v_{i, j}}{h} = 0

Bilinear Interpolation

Use bilinear interpolation to interpolate physical quantities $i\leftarrow [x]\quad j\leftarrow [y]$

f (x) \leftarrow (i + 1 - x) (j + 1 - y) f_{i, j} + (x - i) (j + 1 - y) f_{i + 1, j} + (i + 1 - x) (y - j) f_{i, j + 1} + (x - i) (y - j) f_{i + 1, j + 1}

staggered velocities $x\leftarrow x-0.5$ )

u (x) \leftarrow (i + 1 - x) (j + 1 - y) u_{i, j} + (x - i) (j + 1 - y) u_{i + 1, j} + (i + 1 - x) (y - j) u_{i, j + 1} + (x - i) (y - j) u_{i + 1, j + 1}

Incompressible, Viscous Navier Stokes’ Equations

Equation Fomulation

\begin{matrix} Incompressibility & \nabla \cdot u = 0 \\ Momentum & \frac{\partial u}{\partial t} = g - (u \cdot \nabla) u + μ Δ u - \nabla p \end{matrix}

$\vb{g}$ $(\vb{u}\cdot \grad)\vb{u}$ $\mu\Delta \vb{u}$ $\grad p$ - Internal pressure)

Method of Characteristics: solving a long PDE in steps

$\vb{u}$ $\partial \vb{u}/\partial t = \vb{g}$
$\vb{u}$ $\partial \vb{u}/\partial t = -(\vb{u}\cdot \grad)\vb{u}$
$\vb{u}$ $\partial \vb{u}/\partial t = v\Delta \vb{u}$
$\vb{u}$ $\partial \vb{u}/\partial t = -\grad p$

Step 1: External Acceleration

$u$ $v$

$v$ for verticle velocities <- affected by gravity force)

v_{i, j}^{new} \leftarrow v_{i, j} + Δ t g

Step 2: Advection

(u \cdot \nabla) u = u \cdot \frac{\partial u}{\partial x} + v \cdot \frac{\partial v}{\partial y}

Solving in an Eulerian way can be unstable => use the real meaning of advection: carry physical quantities by velocity (this problem does not exist in Lagrangian methods) => Semi-Lagrangian Method

Semi-Lagrangian Method

Trace a virtual particle backward over time

$\vb{v}$ (verticle))

$\vb{x}_0 \leftarrow (i-0.5, j)$
$\vb{u(x_0)}$
$\vb{x}_1\leftarrow \vb{x_0} -\Delta t\ \vb{u(x_0)}$
$\vb{u(x_1)}$
$u^{\text{new}}_{i,j}\leftarrow u(\vb{x}_1)$

Note: if the velocities are staggered, we need to do staggered bilinear interpolation

=> Subdivide the time step for better tracing

$\Delta t \rightarrow \Delta t/4$ => do 4 times of updates in an original timestep)

$\vb{x}_0\leftarrow(i-0.5, j)$
$\vb{u(x}_0)$
$\vb{x}_1←\vb{x}_0-\frac{1}{3} \Delta t\ \vb u(\vb{x}_0)$
$\vb{u(x}_1)$
$\vb{x}_2\leftarrow \vb x_1- \frac{1}{2} \Delta t\ \vb{u}(\vb{x}_1)$
$\vb{u(x}_2)$
$\vb x_3\leftarrow \vb{x}_2-\frac{1}{2} \Delta t\ \vb{u(x}_2)$
$\vb{u(x}_3)$
$u_{i,j}^{\text{new}}\leftarrow u(\vb{x}_3)$

Step 3: Diffusion

$u$ $v$ )

\begin{aligned} u_{i, j}^{new} \leftarrow u_{i, j} + v Δ t \frac{u_{i - 1, j} + u_{i + 1, j} + u_{i, j - 1} + u_{i, j + 1} - 4 u_{i, j}}{h^{2}} \\ v_{i, j}^{new} \leftarrow v_{i, j} + v Δ t \frac{v_{i - 1, j} + v_{i + 1, j} + v_{i, j - 1} + v_{i, j + 1} - 4 v_{i, j}}{h^{2}} \end{aligned}

$\nu \Delta t$ is large: can be unstable => use smaller sub-steps / implicit integration

Step 4: Pressure Projection

\begin{aligned} u_{i, j}^{n e w} \leftarrow u_{i, j} - \frac{Δ t}{h} (p_{i, j} - p_{i - 1, j}) \\ v_{i, j}^{n e w} \leftarrow v_{i, j} - \frac{Δ t}{h} (p_{i, j} - p_{i, j - 1}) \end{aligned}

The pressure is caused by incompressibilitydivergence free condition $\div\ \vb{u}^{\text{new} }= 0$

\begin{matrix} u_{i + 1, j} + v_{i, j + 1} - u_{i, j} - v_{i, j} = 0 \\ \Rightarrow u_{i + 1, j} - \frac{Δ t (p_{i + 1, j} - p_{i, j})}{h} + v_{i, j + 1} - \frac{Δ t (p_{i, j + 1} - p_{i, j})}{h} - u_{i, j} - \frac{Δ t (p_{i, j} - p_{i - 1, j})}{h} - v_{i, j} - \frac{Δ t (p_{i, j} - p_{i, j - 1})}{h} = 0 \end{matrix}

Eventually got a Poisson equation:

4 p_{i . j} - p_{i - 1, j} - p_{i + 1, j} - p_{i, j - 1} - p_{i, j + 1} = h (- u_{i + 1, j} - v_{i, j + 1} + u_{i, j} + v_{i, j})

With bc:

$p_{i-1,j} = P$
$p_{i-1,j} = p_{i,j}$

$p$ $\vb{u}$ and done

Air and Smoke

Air Simulation

Steps: (can be used in underwater as well)

flow $\vb{u}$
Step 2: Use semi-Lagrangian advect all of the other physical quantites (density / temperature / …)

Typically use Dirichlet boundaries for an open space and Neumann for a container

Water Simulation

Water simulation also involves air-water interface. Water doesn’t occupy the whole space => Need to solve the air-water boundary

Representations:

Volume-of-fluid (VOF) (every grid has a percentage of fluid and air) => not accurate, not tells the interface
A signed distance function defined over grid

Advect: (needs corrections)

Semi-Lagrangian (volume loss)
Level set method (spec. to signed distance function update) (volume loss)

Volume loss affects significantly when perturbation is large